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ABSTRACT 


y  This  report  forms  the  user’s  guide  for  Version  4.0  of  NPSOL,  a  set  of  Fortran  subroutines 
designed  to  minimize  a  smooth  function  subject  to  constraints,  which  may  include  simple  bounds 
on  the  variables,  linear  constraints  and  smooth  nonlinear  constraints.  (NPSOL  may  also  be  used  for 
unconstrained,  bound-constrained  and  linearly  constrained  optimization.)  The  user  must  provide 
subroutines  that  define  the  objective  and  constraint  functions  and  (optionally)  their  gradients.  All 
matrices  are  treated  as  dense,  and  hence  NPSOL  is  not  intended  for  large  sparse  problems. 

NPSOL  uses  a  sequential  quadratic  programming  (SQP)  algorithm,  in  which  the  search  direc¬ 
tion  is  the  solution  of  a  quadratic  programming  (QP)  subproblem.  The  algorithm  treats  bounds, 
iiuear  constraints  and  nonlinear  constraints  separately.  The  Hessian  of  each  QP  subproblcm  is 
a  positive-definite  quasi-Newton  approximation  to  the  Hessian  of  the  Lagrangian  function.  The 
steplength  at  each  iteration  is  required  to  produce  a  sufficient  decrease  in  an  augmented  Lagrangian 
merit  function.  Each  QP  subproblcm  is  solved  using  a  quadratic  programming  pykage  with  several 
features  th^t  improve  the  efficiency  of  an  SQP  algoritl^m. 
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1.  PURPOSE 

NPSOL  is  a  collection  of  Fortran  77  subroutines  designed  to  solve  the  nonlinear  programming 
problem:  the  minimization  of  a  smooth  nonlinear  function  subject  to  a  set  of  constraints  on  the 
variables.  The  problem  is  assumed  to  be  stated  in  the  following  form: 

NP  minimize 

«e»“ 

subject  to 

where  F(x)  (the  objective  function)  is  a  nonlinear  function,  A^,  is  an  x  n  constant  matrix  of 
general  constraints,  and  c(a:)  is  an  m,v-vector  of  nonlinear  constraint  functions.  (The  matrix 
and  the  vector  0(1)  may  be  empty.)  The  objective  function  F  and  the  constraint  functions  are 
eissumed  to  be  smooth,  i.e.,  at  least  twice- continuously  differentiable.  (The  method  of  NPSOL  will 
usually  solve  NP  if  there  are  only  isolated  discontinuities  away  from  the  solution). 

Note  that  upper  and  lower  bounds  are  specified  for  all  the  variables  and  for  all  the  constraints. 
This  form  allows  full  generality  in  specifying  other  types  of  constraints.  In  particular,  the  i-th 
constraint  may  be  defined  as  an  equality  by  setting  /j  =  If  certain  bounds  are  not  present,  the 
associated  elements  of  /  or  u  can  be  set  to  special  values  that  will  be  treated  as  —00  or  -foe. 

If  there  are  no  nonlinear  constraints  in  NP  and  F  is  linear  or  quadratic,  the  QPSOL  or  LSSOL 
packages  (Gill  et  al.,  1984a,  1986a)  will  generally  be  more  efficient  than  NPSOL.  If  the  problem  is 
large  and  sparse,  the  MINOS  package  (Murtagh  and  Saunders,  1982,  1983)  should  be  used,  since 
NPSOL  treats  all  matrices  as  dense. 

The  user  must  supply  an  initial  estimate  of  the  solution  to  NP,  and  subroutines  that  define 
F{x),  c(i),  and  as  many  first  partial  derivatives  as  possible;  unspecified  derivatives  are  approxi¬ 
mated  by  finite-differences. 

NPSOL  is  based  on  subroutines  from  Version  1.0  of  the  LSSOL  constrained  linear  least-squares 
package;  the  documentation  of  LSSOL  (Gill  et  al.,  1986a)  should  be  consulted  in  conjunction  with 
this  report.  NPSOL  contains  approximately  15,000  lines  of  ANSI  (1977)  Standard  Fortran,  of  which 
47%  are  comments. 
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2.  DESCRIPTION  OF  THE  ALGORITHM 

Here  we  briefly  summarize  the  main  features  of  the  method  of  NPSOL.  Where  possible,  explicit 
reference  is  made  to  the  names  of  variables  that  are  parameters  of  subroutine  NPSOL  or  appear  in 
the  printed  output. 


At  a  solution  of  NP,  some  of  the  constraints  will  be  active,  i.e.,  satisfled  exactly.  An  active 
simple  bound  constraint  implies  that  the  corresponding  variable  is  fixed  at  its  bound,  and  hence 
the  variables  are  partitioned  into  fixed  and  free  variables.  Let  C  denote  the  m  x  n  matrix  of 
gradients  of  the  active  general  linear  and  nonlinear  constraints.  The  number  of  fixed  variables  will 
be  denoted  by  n^x,  with  (npR  =  n  —  npx)  the  number  of  free  variables.  The  subscripts  “FX” 
and  “FR”  on  a  vector  or  matrix  will  denote  the  vector  or  matrix  composed  of  the  components 
corresponding  to  fixed  or  free  variables. 

A  point  I  is  a  first-order  Kuhn-Tucker  point  for  NP  (see,  e.g.,  Powell,  1974)  if  the  following 
conditions  hold: 


(i)  X  is  feasible; 

(ii)  there  exist  vectors  ^  and  A  (the  Lagrange  multiplier  vectors  for  the  bound  and  general 
constraints)  such  that 

g  =  C^\  +  C,  (1) 

where  g  is  the  gradient  of  F  evaluated  at  x,  and  =  0  if  the  j-th  variable  is  free. 

(iii)  The  Lagrange  multiplier  corresponding  to  an  inequality  constraint  active  at  its  lower 
bound  must  be  non-negative,  and  non-positive  for  an  inequality  constraint  active  at 
its  upper  bound. 

Let  Z  denote  a  matrix  whose  columns  form  a  basis  for  the  set  of  vectors  orthogonal  to  the 
rows  of  Cfr;  i.e.,  Cfn,Z  =  0.  An  equivalent  statement  of  the  condition  (1)  in  terms  of  Z  is 


=  0. 


The  vector  Z^g^n  is  termed  the  projected  gradient  of  F  at  x.  Certain  additional  conditions  must 
be  satisfied  in  order  for  a  first-order  Kuhn-Tucker  point  to  be  a  solution  of  NP  (see,  e.g.,  Powell, 
1974). 

The  method  of  NPSOL  4.0  is  a  sequential  quadratic  programming  (SQP)  method.  For  an 
overview  of  SQP  methods,  see,  for  example,  Fletcher  (1981),  Gill,  Murray  and  Wright  (1981)  and 
Powell  (1983). 

The  basic  structure  of  NPSOL  involves  major  and  minor  iterations.  The  major  iterations 
generate  a  sequence  of  iterates  {ifc}  that  converge  to  a:*  a  first-order  Kuhn-Tucker  point  of  NP.  At 
a  typical  major  iteration,  the  new  iterate  x  is  defined  by 


X  =  X  +  ap, 


(2) 


where  x  is  the  current  iterate,  the  non-negative  scalar  a  is  the  step  length,  and  p  is  the  search 
direction.  (For  simplicity,  we  shall  always  consider  a  typiced  iteration  and  avoid  reference  to  the 
index  of  the  iteration.)  Also  associated  with  each  major  iteration  are  estimates  of  the  Lagrange 
multipliers  and  a  prediction  of  the  active  set. 

The  search  direction  p  in  (2)  is  the  solution  of  a  quadratic  programming  subproblera  of  the 

form 

minimize  g^p  + 

»  ^ 


subject  to 


(3) 
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where  g  is  the  gradient  of  F  at  x,  the  matrix  H  is  &  positive-definite  quasi-Newton  api)roximation 
to  the  Hessian  of  the  Lagrangian  function  (see  Section  2.3),  and  i4jv  is  the  Jacobian  matrix  of  c 
evaluated  at  x.  (Finite-difference  estimates  may  be  used  for  g  and  >1;^;  see  the  optional  parameter 
“Derivative  Level”  in  Section  5.2.)  Let  t  in  NP  be  partitioned  into  three  sections;  tg,  II  and 
/jv,  corresponding  to  the  bound,  linear  and  nonlinear  constraints.  The  vector  I  in  (3)  is  similarly 
partitioned,  and  is  defined  as 

Ig  “  Ig  “  X,  j4^x,  and 

where  c  is  the  vector  of  nonlinear  constraints  evaluated  at  x.  The  vector  u  is  defined  in  an  analogous 
fashion. 

The  estimated  Lagrange  multipUers  at  each  major  iteration  are  the  Lagrange  multipliers  from 
the  subproblem  (3)  (and  similarly  for  the  predicted  active  set).  (The  numbers  of  bounds,  general 
linear  and  nonlinear  constraints  in  the  QP  active  set  are  the  quantities  “End”,  “Lin”  and  “Nln” 
in  the  printed  output  of  NPSOL.)  In  NPSOL,  (3)  is  solved  using  subroutines  from  Version  1.0  of 
the  LSSOL  package  (Gill  et  ai.,  1986a).  Since  solving  a  quadratic  program  is  itself  an  iterative 
procedure,  the  minor  iterations  of  NPSOL  are  the  iterations  of  LSSOL.  (More  details  about  solving 
the  subproblem  are  given  in  Section  2.1.) 

Certain  matrices  associated  with  the  QP  subproblem  are  relevant  in  the  major  iterations.  Let 
the  subscripts  “FX”  and  “FR”  refer  to  the  predicted  fixed  and  free  variables,  and  let  C  denote  the 
m  X  n  matrix  of  gradients  of  the  general  linear  and  nonlinear  constraints  in  the  predicted  active 
set.  First,  we  have  available  the  TQ  factorization  of  Crg'- 

C,gQ,g  =  (  0  T  ),  (4) 

where  T  is  a  nonsingular  m  x  m  reverse- triangular  matrix  (i.e.,  tij  =  0  if  t  -f  jf  <  m),  and  the 
non-singular  ripR  x  matrix  is  the  product  of  orthogonal  transformations  (see  Gill  et  aJ., 
1984a).  Second,  we  have  the  upper-triangular  Cholesky  factor  R  of  the  transformed  and  re-ordered 
Hessian  matrix 

R^R  =  Hq  =  Q^HQ,  (5) 

where  H  is  the  Hessian  H  v.rith  rows  and  columns  permuted  so  that  the  free  variables  are  first,  smd 
Q  is  the  n  X  n  matrix 


with  /px  the  identity  matrix  of  order  Upx-  If  the  columns  of  Qrg  are  partitioned  so  that 

Qrg={Z  Y), 


the  Tip  (rip  =  Tipp  —  m)  columns  of  Z  form  a  basis  for  the  null  space  of  CpR.  The  matrix  Z  is  used 
to  compute  the  projected  gradient  Z^g^g  at  the  current  iterate.  (The  values  “Nz”,  “Norm  Gf”,  and 
“Norm  Gz”  printed  by  NPSOL  give  Up  and  the  norms  of  gfg  and  Z^g^g.) 

A  theoretical  characteristic  of  SQP  methods  is  that  the  predicted  active  set  from  the  QP 
subproblem  (3)  is  identical  to  the  correct  active  set  in  a  neighborhood  of  x*  In  NPSOL,  this  feature 
is  exploited  by  using  the  QP  active  set  from  the  previous  iteration  as  a  prediction  of  the  active 
set  for  the  next  QP  subproblem,  which  leads  in  practice  to  optimality  of  the  subproblems  in  only 
one  iteration  as  the  solution  is  approached.  Separate  treatment  of  bound  and  linear  constraints  in 
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Once  p  lias  been  computed,  the  major  iteration  proceeds  by  determining  a  stcplength  a  that 
produces  a  “sufficient  decrease”  in  an  augmented  Lagrangian  merit  function  (see  Section  2.2). 
Finally,  the  approximation  to  the  transformed  Hessian  matrix  JIq  is  updated  using  a  modified 
BFGS  quasi-Newton  update  (see  Section  2.3)  to  incorporate  new  curvature  information  obtained 
in  the  move  from  x  to  x. 

On  entry  to  NPSOL,  an  iterative  procedure  from  the  LSSOL  package  is  executed,  starting  with 
the  user-provided  initial  point,  to  find  a  point  that  is  feasible  with  respect  to  the  bounds  and  linear 
constraints  (using  the  tolerance  specified  by  “Linear  Feasibility  Tolerance”;  see  Section  5.2). 
If  no  feasible  point  exists  for  the  bound  and  linear  constraints,  NP  has  no  solution  and  NPSOL 
terminates.  Otherwise,  the  problem  functions  will  thereafter  be  evaluated  only  at  points  that  are 
feasible  with  respect  to  the  bounds  and  linear  constraints.  The  only  exception  involves  variables 
whose  bounds  differ  by  an  amount  comparable  to  the  finite-difference  interval  (see  the  discussion 
of  “Difference  Interval”  in  Section  5.2).  In  contrast  to  the  bounds  and  linear  constraints,  it 
must  be  emphasized  that  the  nonlinear  constraints  will  not  generally  be  satisfied  until  an  optimal 
point  is  reached. 

Facilities  are  provided  to  check  whether  the  user-provided  gradients  appear  to  be  correct  (see 
the  optional  parameter  “Verify”  in  Section  5.2).  In  general,  the  check  is  provhled  at  the  first 
point  that  is  feasible  with  respect  to  the  linear  constraints  and  bounds.  However,  the  user  may 
request  that  the  check  be  performed  at  the  initial  point. 

In  summary,  the  method  of  NPSOL  first  determines  a  point  that  satisfies  the  bound  and  linear 
constraints.  Thereafter,  each  iteration  includes:  (a)  the  solution  of  a  quadratic  programming 
subproblem;  (b)  a  linesearch  with  an  augmented  Lagrangian  merit  function;  and  (c)  a  quasi- 
Newton  update  of  the  approximate  Hessian  of  the  Lagrangian  function.  These  three  procedures 
are  described  in  more  detail  in  the  next  three  subsections. 

2.1.  Solution  of  the  quadratic  programming  subproblem 

The  search  direction  p  is  obtained  by  solving  (3)  using  subroutines  from  the  LSSOL  package  (Gill 
ef  at.,  1986a),  which  was  specifically  designed  to  be  used  within  an  SQP  algorithm  for  nonlinear 
programming. 

The  method  of  LSSOL  is  a  two-phase  (primal)  quadratic  programming  method.  The  two 
phases  of  the  method  are:  finding  an  initial  feasible  point  by  minimizing  the  sum  of  infeasibilities 
(the  feasibility  phase),  and  minimizing  the  quadratic  objective  function  within  the  feasible  region 
(the  optimality  phase).  The  computations  in  both  phases  are  performed  by  the  same  subroutines. 
The  two-phase  nature  of  the  algorithm  is  reflected  by  changing  the  function  being  minimized  from 
the  sum  of  infeasibilities  to  the  quadratic  objective  function. 

In  general,  a  quadratic  program  must  be  solved  by  iteration.  Let  p  denote  the  current  estimate 
of  the  solution  of  (3);  the  new  iterate  p  is  defined  by 

p  =  p  -f  <rd,  (7) 

where,  as  in  (2),  <7  is  a  non-negative  step  length  and  d  is  a  search  direction. 

At  the  beginning  of  each  iteration  of  LSSOL,  a  working  set  is  defined  of  constraints  (general 
and  bound)  that  are  satisfied  exactly.  The  vector  d  is  then  constructed  so  that  the  values  of 
constraints  in  the  working  set  remain  unaltered  for  any  move  along  d.  For  a  bound  constraint  in 
the  working  set,  this  property  is  achieved  by  setting  the  corresponding  component  of  d  to  zero, 
i.e.,  by  fixing  the  variable  at  its  bound.  As  before,  the  subscripts  “FX”  and  “FR”  denote  selection 
of  the  components  associated  with  the  fixed  and  free  variables. 
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Let  C  denote  the  submatrix  of  rows  of 


corresponding  to  general  constraints  in  the  working  set.  The  general  constraints  in  the  working 
set  will  remain  unaltered  if 

Cfr^ph  =  0»  (8) 

which  is  equivalent  to  defining  dpR  as 

dpB  =  Zdg  (9) 

for  some  vector  d2,  where  Z  is  the  matrix  associated  with  the  TQ  factorization  (4)  of  CpR. 

The  definition  oldz  in  (9)  depends  on  whether  the  current  p  is  feasible.  If  not,  d^  is  zero  except 
for  a  component  7  in  the  j-th  position,  where  j  and  7  are  chosen  so  that  the  sum  of  infeasibilities 
is  decreasing  along  d.  (For  further  details,  see  Gill  et  aJ.,  1986a.)  In  the  feasible  case,  d^  satisfies 
the  equations 

R'^2R^d2  =  -Z%^,  (10) 

where  is  the  Cholesky  factor  of  Z  and  q  is  the  gradient  of  the  quadratic  objective  function 

(9  =  3  +  Rp)-  (The  vector  is  the  projected  gradient  of  the  QP.)  With  (10),  p  +  d  is  the 

miniraizer  of  the  quadratic  objective  function  subject  to  treating  the  constraints  in  the  working 
set  as  equalities. 

If  the  QP  projected  gradient  is  zero,  the  current  point  is  a  constrained  stationary  point  in 
the  subspace  defined  by  the  working  set.  During  the  feasibility  phase,  the  projected  gradient  will 
usually  be  zero  only  at  a  vertex  (although  it  may  vanish  at  non-vertices  in  the  presence  of  constraint 
dependencies).  During  the  optimality  phase,  a  zero  projected  gradient  implies  that  p  minimizes 
the  quadratic  objective  function  when  the  constraints  in  the  working  set  are  treated  as  equalities. 
In  either  case,  Lagrange  multipliers  are  computed.  Given  a  positive  constant  6  of  the  order  of 
the  machine  precision,  the  Lagrange  multiplier  fij  corresponding  to  an  inequality  constraint  in  the 
working  set  at  its  upper  bound  is  said  to  be  optimal  if  pj  <  6  when  the  j-th  constraint  is  at  its 
upper  bound,  or  if  py  >  —6  when  the  associated  constraint  is  at  its  lower  bound.  If  any  multiplier  is 
non-optimal,  the  current  objective  function  (either  the  true  objective  or  the  sum  of  infeasibilities) 
can  be  reduced  by  deleting  the  corresponding  constraint  from  the  working  set. 

If  optimal  multipliers  occur  during  the  feasibility  phase  and  the  sum  of  infeasibilities  is  nonzero, 
no  feasible  point  exists.  The  QP  algorithm  will  then  continue  iterating  to  determine  the  minimum 
sum  of  infeasibilities.  At  this  point,  the  Lagrange  multiplier  pj  will  satisfy  —(1  +  6)  <  pj  <  S  for 
an  inequality  constraint  at  its  upper  bound,  and  -5  <  <  1  -I-  5  for  an  inequality  at  its  lower 

bound.  The  Lagrange  multiplier  for  an  equality  constraint  will  satisfy  |pj|  <  1  -t-  6. 

The  choice  of  step  length  a  in  the  QP  iteration  (7)  is  based  on  remaining  feasible  with  respect 
to  the  satisfied  constraints.  During  the  optimality  phase,  if  p  -f  d  is  feasible,  a  will  be  taken  as 
unity.  (In  this  case,  the  projected  gradient  at  p  will  be  zero.)  Otherwise,  a  is  set  to  the  step 
to  the  “nearest”  constraint,  which  is  added  to  the  working  set  at  the  next  iteration. 

Each  change  in  the  working  set  leads  to  a  simple  change  to  Gfr:  if  the  status  of  a  general 
constraint  changes,  a  row  of  Ctr  is  altered;  if  a  bound  constraint  enters  or  leaves  the  working  set, 
a  column  of  Crr  changes.  Explicit  representations  are  recurred  of  the  matrices  T,  and  i?,  and 
of  the  vectors  Q^q  and  Q^g. 
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2.2.  The  merit  function 

After  computing  the  search  direction  as  described  in  Section  2.1,  each  major  iteration  proceeds  by 
determining  a  steplength  a  in  (2)  that  produces  a  “sufficient  decrease”  in  the  augmented  Lagrangian 
merit  function 

£(i,A,.s)  =  F{x)  -  Ai(c<(i)  -  Si)  +  ^  -  Si)^  (11) 

i  i 

where  z,  A  and  s  vary  during  the  linesearch.  The  summation  terms  in  (11)  involve  only  the 
nonlinear  constraints.  The  vector  A  is  an  estimate  of  the  Lagrange  multipliers  for  the  nonlinear 
constraints  of  NP.  The  non-negative  slack  variables  {s^}  allow  nonlinear  inequality  constraints  to 
be  treated  without  introducing  discontinuities.  The  solution  of  the  QP  subproblem  (3)  provides  a 
vector  triple  that  serves  as  a  direction  of  search  for  the  three  sets  of  variables.  The  non-negative 
vector  p  of  penalty  parameters  is  initialized  to  zero  at  the  beginning  of  the  first  major  iteration. 
Thereafter,  selected  components  are  increased  whenever  necessary  to  ensure  descent  for  the  merit 
function.  Thus,  the  sequence  of  norms  of  p  (the  printed  quantity  “Penalty”;  see  Section  6)  is 
generally  non-decreasing,  although  each  pi  may  be  reduced  a  limited  number  of  times. 

The  merit  function  (11)  and  its  global  convergence  properties  are  deseribed  in  Gill  et  al. 
(1986b). 


2.3.  The  quasi-Newton  update 

The  matrix  H  in  (3)  is  a  positive-definite  quasi-Newton  approximation  to  the  Hessian  of  the  La¬ 
grangian  function.  (For  a  review  of  quasi-Newton  methods,  see  Dennis  and  Schnabel,  1983.)  At  the 
end  of  each  major  iteration,  a  new  Hessian  approximation  ff  is  defined  as  a  rank-two  modification 
of  ff.  In  NPSOL,  the  BFGS  quasi-Newton  update  is  used: 

where  s  =  x  —  z  (the  change  in  z). 

In  NPSOL,  ff  is  required  to  be  positive  definite.  If  ff  is  positive  definite,  ff  as  defined  by  (12) 
will  be  positive  definite  if  and  only  if  y^s  is  positive  (see,  e.g.,  Dennis  and  More,  1977).  Ideally,  y 
in  (12)  would  be  taken  as  j/t,  the  change  in  gradient  of  the  Lagrangian  function 

yt  =  5  -  -  y  +  (13) 


where  denotes  the  QP  multipliers  asociated  with  the  nonlinear  constraints  of  the  original 
problem.  If  y^s  is  not  sufficiently  positive,  an  attempt  is  made  to  perform  the  update  with  a  vector 
y  of  the  form 

mu 

y  =  Vt  +  ]^Wi(o<(i)ci(x)  -  a<(z)ci(z)), 

»=i 

where  Wj  >  0.  If  no  such  vector  can  be  found,  the  update  is  performed  with  a  scaled  j/i,;  in  this 
case,  “M”  is  printed  to  indicate  that  the  update  was  modified. 

Rather  than  modifying  ff  itself,  the  Cholcsky  factor  of  the  transformed  Hessian  Hq  (4)  is 
updated,  where  Q  is  the  matrix  from  (3)  associated  with  the  active  set  of  the  QP  subproblem.  The 
update  (12)  is  equivalent  to  the  following  update  to  Hq'. 


Hq  —  ffq  -  _  HqSqS^Hq  -f-  „  yoyq, 

where  =  Q^y,  and  Sq  =  Q^s.  This  update  may  be  expressed  as  a  rank-one  update  to 
Dennis  and  Schnabel,  1981). 

Full  details  concerning  the  Hessian  update  are  given  in  Gill  et  al.  (1986c). 


(14) 
ff  (see 
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3.  SPECIFICATION  OF  SUBROUTINE  NPSOL 

The  formal  specification  of  NPSOL  is  the  following: 


SUBROUTINE  NPSOL  (  N.  NCLIN,  NCNLN,  NROWA,  NROWJ,  NROWR, 
A,  BL.  BU, 

CONFUN,  OBJFUN, 

INFORM,  ITER,  ISTATE, 

C,  CJAC,  CLAMDA,  OBJF.  GRAD,  R,  X, 

IW,  LENIW.  W,  LENW  ) 


INTEGER 

INTEGER 

REAL 

REAL 


EXTERNAL 


N,  NCLIN,  NCNLN, 

NROWA,  NROWJ,  NROWR,  INFORM,  ITER,  LENIW,  LENW 
ISTATE(N+NCLIN+NCNLN) ,  IW(LENIW) 

OBJF 

A(NROWA,*),  BL(N+NCLIN+NCNLN) ,  BU(N+NCLIN+NCNLN) , 
C(*),  CJAC (NROWJ,*),  CLAMDA (N+NCLIN+NCNLN),  GRAD(N), 
RCNROWR,*),  X(N),  W(LENW) 

CONFUN,  OBJFUN 


Note:  Here  and  elsewhere,  the  specification  of  a  parameter  as  REAL  should  be  interpreted  as  working 
precision,  which  may  be  DOUBLE  PRECISION  in  some  installations. 


3.1.  Formal  parameters 

N  (Input)  The  number  of  variables  in  the  problem,  i.e.,  the  dimension  of  X.  (N  must 

be  positive.) 

NCLIN  (Input)  The  number  of  general  linear  constraints  in  the  problem.  '{NCLIN  may  be 

zero.) 

NCNLN  (Input)  The  number  of  nonlinear  constraints  in  the  problem.  (NCNLN  may  be  zero.) 

NROWA  (Input)  The  declared  row  dimension  of  the  array  A.  NROWA  must  be  at  least  1  and 

at  least  NCLIN. 

NROWJ  (Input)  The  declared  row  dimension  of  the  array  CJAC.  NROWJ  must  be  at  least  1 

and  at  least  NCNLN. 

NROWR  (Input)  The  declared  row  dimension  of  the  array  R.  NROWR  must  be  at  least  N. 

A  (Input)  A  real  array  of  declared  dimension  (NROWA,*),  where  the  second  dimension 

must  be  at  least  N.  A  contains  the  matrix  Aj,  of  general  linear  constraints  in  the 
problem  specification  NP  (Section  1).  The  t-th  row  of  A,  i  =  1  to  NCLIN,  contains  the 
coefficients  of  the  t-th  general  linear  constraint.  If  NCLIN  is  zero,  A  is  not  accessed 
and  may  be  dimensioned  (1,1). 

BL  (Input)  A  real  array  of  dimension  at  least  N -|- NCLIN -|- NCNLN  that  contains  the  lower 

bounds  for  all  the  constraints,  in  the  following  order  (which  is  also  observed  for  BU, 
CLAMDA  and  ISTATE).  The  first  N  elements  of  BL  contain  the  lower  bounds  on  the 
variables.  If  NCLIN  >  0,  the  next  NCLIN  elements  of  BL  contain  the  lower  bounds  for 
the  general  hnear  constraints.  If  NCNLN  >  0,  the  next  NCNLN  elements  of  BL  contain 
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BU 


CONFUN 


OBJFUN 


INFORM 


<  0  The  user  has  set  MODE  to  this  negative  value  in  CONFUN  or  OBJFUN  (see 
Section  4). 

0  The  iterates  have  converged  to  a  point  X  that  satisfies  the  first-order 
Kuhn-Tucker  conditions  to  the  accuracy  requested  by  the  optional  pa¬ 
rameter  Optimality  Tolerance  (see  Section  5.2),  i.e.,  the  projected  grac 
diont  and  activ*'  constraint  re.siduals  are  negligible  at  X. 

1  The  final  iterate  X  satisfies  the  first-order  Kuhn-Tucker  conditions  to  the 
accuracy  requested,  but  the  sequence  of  iterates  has  not  yet  converged. 
NPSOL  was  terminated  because  no  further  improvement  could  be  made 
in  the  merit  function. 

2  No  feasible  point  could  be  found  for  the  linear  constraints  and  bounds. 
The  problem  has  no  feasible  solution.  See  Section  7  for  further  com¬ 
ments. 

3  No  feasible  point  could  be  found  for  the  nonlinear  constraints.  The  prob¬ 
lem  may  have  no  feasible  solution.  See  Section  7  for  further  comments. 

4  The  limiting  number  of  iterations  (determined  by  the  optional  parameter 
Major  Iteration  Limit;  sec  Section  5.2)  has  been  reached. 

6  X  docs  not  satisfy  the  first-order  Kuhn-Tucker  conditions,  and  no  im¬ 
proved  point  for  the  merit  function  could  be  found  during  the  final  line 
search. 

7  The  user-provided  derivatives  of  the  objective  function  and/or  nonbnear 
constraints  appear  to  be  incorrect. 

9  An  input  parameter  is  invalid. 

(Output)  The  number  of  major  iterations  performed. 


the  lower  bounds  for  the  nonlinear  constraints.  In  order  for  the  problem  specification 
to  be  meaningful,  it  is  re<piired  that  BL(ji)  <  BU(j)  for  all  j.  To  specify  a  non-existent 
lower  bound  (i.e.,  tj  =  — oc),  the  value  used  must  satisfy  BL(j)  <  -BIGBND,  where 
BIGBND  is  the  vahie  of  the  optional  parameter  Infinite  Bound,  whose  default  value 
is  10*®  (see  Section  5.2).  To  specify  the  y-th  constraint  as  an  equality,  the  user  must 
set  BL(j)  =  BU(j)  /3,  say,  where  |/3|  <  BIGBND. 

(Input)  A  real  array  of  dimension  at  least  N -f  NCLIN-f  NCNLN  that  contains  the  upper 
bounds  for  all  the  constraints,  in  the  same  order  described  above  for  BL.  To  specify  a 
non-existent  ui)por  bound  (i.e.,  Uj  =  oo),  the  value  used  must  satisfy  BU(y  )  >  BIGBND. 

(User- defined  subroutine)  The  name  of  a  subroutine  that  calculates  the  vector 
c(x)  of  nonlinear  constraint  functions  and  (optionally)  its  Jacobian  for  a  specified 
n-vector  x.  CONFUN  must  be  declared  as  EXTERNAL  in  the  routine  that  calls  NPSOL. 
For  a  detailed  description  of  CONFUN,  sec  Section  4.2. 

(User-defined  subroutine)  The  name  of  a  subroutine  that  calculates  the  objective 
function  F{x)  and  (optionally)  its  gradient  for  a  specified  n-vector  x.  OBJFUN  must 
be  declared  as  EXTERNAL  in  the  routine  that  calls  NPSOL.  For  a  detailed  description 
of  OBJFUN,  see  Section  4.1. 

(Output)  An  integer  that  indicates  the  result  of  NPSOL.  (A  short  description  of 
INFORM  is  printed  if  Major  Print  Level  >  0.)  The  possible  values  of  INFORM  are: 

INFORM  Meaning 


i 


e 


u 


C, 

A 


i.) 


ITER 
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ISTATE 


(Input)  An  integer  array  of  dimension  at  least  N  +  NCLIN  +  NCNLN.  ISTATE  need  not 
be  initialized  if  NPSOL  is  called  with  a  Cold  Start  (the  default  option;  sec  Section 
5.2).  The  ordering  of  ISTATE  is  the  same  as  that  described  above  for  BL,  i.c.,  the 
first  N  components  of  ISTATE  refer  to  the  upper  and  lower  bounds  on  the  variables, 
components  N  +  1  through  N  +  NCLIN  refer  to  the  upper  and  lower  bounds  on 
and  components  N  +  NCLIN  +  1  through  N  +  NCLIN  +  NCNLN  refer  to  the  upper  and 
lower  bounds  on  c(a;).  When  a  Warm  Start  option  is  chosen,  the  components  of 
ISTATE  corresponding  to  the  bounds  and  linear  constraints  define  the  initial  working 
set  for  the  procedure  that  finds  a  feasible  point  for  the  linear  constraints  and  bounds. 
The  active  set  at  the  conclusion  of  this  procedure  and  the  components  of  ISTATE 
corresponding  to  nonbnear  constraints  then  define  the  initial  working  set  for  the  first 
QP  subproblem.  Possible  values  for  ISTATE(j)  are 


ISTATE(i) 

0 

1 

2 

3 


Meaning 

The  corresponding  constraint  is  not  in  the  initial  QP  working  set. 

This  inequality  constraint  should  be  in  the  working  set  at  its  lower  bound. 
This  inequality  constraint  should  be  in  the  working  set  at  its  upper 
bound. 

This  equality  constraint  should  be  in  the  initial  working  set.  This  value 
must  not  be  specified  unless  BL(j)  =  BU(j).  The  values  1,  2  or  3  all  have 
the  same  effect  when  BL(j)  =  BU(j). 


Other  values  of  ISTATE  are  also  acceptable.  In  particular,  if  NPSOL  has  been  called 
previously  with  the  same  values  of  N,  NCLIN  and  NCNLN,  ISTATE  already  contains  sat¬ 
isfactory  values.  If  necessary,  NPSOL  will  override  the  user’s  specification  of  ISTATE, 
so  that  a  poor  choice  will  not  cause  the  algorithm  to  fail. 


(Output)  If  NPSOL  exits  with  INFORM  =  0  or  1,  the  values  in  the  array  ISTATE 
correspond  to  the  active  set  of  the  final  QP  subproblem,  and  are  a  prcdictioa  of  th^ 
status  of  the  constraints  at  the  solution  of  the  problem.  Otherwise,  ISTATE  indicates 
the  composition  of  the  QP  working  set  at  the  final  iterate.  The  significance  of  each 
possible  value  of  ISTATE(j)  is  as  follows: 

ISTATE(j)  Meaning 

—2  This  constraint  violates  its  lower  bound  by  more  than  the  feasibility 
tolerance  (see  the  optional  parameters  Lineeu:  Feasibility  Tolerance 
and  Nonlinear  Feasibility  Tolerance  in  Section  5.2).  This  value  can 
occur  only  when  no  feasible  point  can  be  found  for  a  QP  subproblem. 

-1  This  constraint  violates  its  upper  bound  by  more  than  the  appropri¬ 
ate  feasibility  tolerance  (see  the  optional  parameters  Linear  Feasi¬ 
bility  Tolerance  and  Nonlinear  Feasibility  Tolerance  in  Section 
5.2).  This  value  can  occur  only  when  no  feasible  point  can  be  found  for 
a  QP  subproblem. 

0  The  constraint  is  satisfied  to  within  the  feasibility  tolerance,  but  is  not 
in  the  working  set. 

1  This  inequality  constraint  is  included  in  the  QP  working  set  at  its  lower 
bound. 

2  This  inequality  constraint  is  included  in  the  QP  working  set  at  its  upper 
bound. 
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3  This  constraint  is  included  in  the  QP  working  set  as  an  equality.  This 
value  of  ISTATE  can  occur  only  when  BL(j)  =  BU(j). 

C  (Output)  A  real  array  of  dimension  at  least  NCNLN.  If  NCNLN  =  0,  C  is  not  accessed, 

and  may  then  be  declared  to  be  of  dimension  (1),  or  the  actual  parameter  may  be 
any  convenient  array.  If  NCNLN  >  0,  C  contains  the  values  of  the  nonlinear  constraint 
functions  Cj,  t  =  1  to  NCNLN,  at  the  final  iterate. 

CJAC  (Input)  A  real  array  of  dimension  (NROWJ,*),  where  the  second  dimension  must  be 

at  least  N.  If  NCNLN  =  0,  CJAC  is  not  accessed,  and  may  then  be  declared  to  be  of 
d  lension  (1,1),  or  the  actual  parameter  may  be  any  convenient  array. 

In  general,  CJAC  need  not  be  initialized  before  the  call  to  NPSOL.  However,  if  Deriva¬ 
tive  Level  =  3,  the  user  may  optionally  set  the  constant  elements  of  CJAC  (see  Sec¬ 
tion  4.3).  Such  constant  elements  need  not  be  re-assigned  on  subsequent  calls  to 
CONFUN. 

(Output)  If  NCNLN  >  0,  CJAC  contains  the  Jacobian  matrix  of  the  nonlinear  con¬ 
straint  functions  at  the  final  iterate,  i.e.,  CJAC(t,_7)  contains  the  partial  derivative  of 
the  t-th  constraint  function  with  respect  to  the  y-th  variable,  i  =  1  to  NCNLN,  j  =  I 
to  N.  (See  the  discussion  of  CJAC  under  CONFUN  in  Section  4.2.) 

CLAMDA  (Input)  A  real  array  of  dimension  at  least  N  -f-  NCLIN  -f  NCNLN.  CLAMDA  need  not  be 
initialized  if  NPSOL  is  called  with  the  (default)  Cold  Start  option.  With  the  Warm 
Start  option,  CLAMDA  must  contain  a  multiplier  estimate  for  each  nonlinear  constraint 
with  a  sign  that  matches  the  status  of  the  constraint  specified  by  the  ISTATE  array 
(as  above).  The  ordering  of  CLAMDA  is  the  same  as  that  given  above  for  BL.  If 
the  j-th  constraint  is  defined  as  “inactive”  by  the  initial  value  of  the  ISTATE  array, 
CLAMDA(y)  should  be  zero;  if  the  j-th  constraint  is  an  ineqxiality  active  at  its  lower 
bound,  CLAMDA(j)  should  be  non-negative;  if  the  j-th  constraint  is  an  iftequality  active 
at  its  upper  bound,  CLAMDA(j)  should  be  non-positive. 

(Output)  CLAMDA  gives  the  QP  multipliers  from  the  last  QP  subproblem.  CLAMDA(j) 
should  be  non-negative  if  ISTATE(j)  =  1  and  non-positive  if  ISTATE(j)  =  2. 

OBJF  (Output)  The  value  of  the  objective  function  F(x)  at  the  final  iterate. 

OBJGRD  (Output)  A  real  array  of  dimension  at  least  N  that  contains  the  objective  gradient 

(or  its  finite-difference  approximation)  at  the  final  iterate. 

R  (Input)  A  real  array  of  declared  dimension  (NROWR,*),  where  the  second  dimension 

must  be  at  least  N.  R  need  not  be  initialized  if  NPSOL  is  called  with  a  Cold  Start 
option  (the  default),  and  will  be  taken  as  the  identity.  With  a  Warm  Start,  R  must 
contain  the  upper-triangular  Cholesky  factor  of  the  initial  approximation  of  the  Hes¬ 
sian  of  the  Lagrangian  function,  with  the  variables  in  the  natural  order.  Elements  not 
in  the  upper-triangular  part  of  R  are  assumed  to  be  zero  and  need  not  be  assigned. 

(Output)  If  Hessian  =  No  (the  default;  see  Section  5.2),  R  contains  the  upper- 
triangular  Cholesky  factor  of  Q'^HQ,  an  estimate  of  the  transformed  and  re-ordered 
Hessian  of  the  Lagrangian  at  X  (see  (5)  in  Section  2).  If  Hessian  »  Yes,  R  contains 
the  upper- triangular  Cholesky  factor  of  H,  the  approximate  (untransformed)  Hessian 
of  the  Lagrangian,  with  the  variables  in  the  natural  order. 
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I  (Input)  A  roal  array  of  dimension  at  least  N.  X  must  contain  an  initial  estimate  of 

the  solution. 

(Output)  X  contains  the  final  estimate  of  the  solution. 


3.2.  Workspace  parameters 

IW  (Input)  An  integer  array  of  dimension  LENIV  that  provides  integer  workspace  for 

NPSOL. 

LENIW  (Input)  The  dimension  of  IW.  LENIW  must  be  at  least  3  N  +  NCLIN  +  2  NCNLN. 

W  (Input)  A  real  array  of  dimension  LENV  that  provides  real  workspace  for  NPSOL. 

LENW  (Input)  The  dimension  of  W.  If  there  are  no  general  linear  constraints  and  no  nonlin¬ 

ear  constraints  (i.e.,  NCLIN  =  0  and  NCNLN  =  0),  LENW  must  be  at  least  20  N.  If  there  are 
no  nonlinear  constraints  (i.e.,  NCNLN  =  0),  LENW  must  be  at  least  2  N^  -|-  20  N  -|- 11  NCLIN. 
Otherwise,  LENW  must  be  at  least  +  N+NCLIN  +  2N+NCNLN  +  20  N+  11  NCLIN  -f 
21  NCNLN. 

If  Major  Print  L«v«l  >  0,  the  required  amounts  of  workspace  are  printed.  As  an  alternative 
to  computing  LENIW  and  LENW  from  the  formulas  given  above,  the  user  may  prefer  to  obtain 
appropriate  values  from  the  output  of  a  preliminary  run  with  a  positive  value  of  Major  Print 
Level  and  LENIW  and  LENW  set  to  1.  (NPSOL  wUl  then  terminate  with  INFORM  =  9.) 
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4.  USER  SUPPLIED  SUBROUTINES 

The  user  must  provide  subroutines  that  define  the  objective  function  and  nonhnear  constraints. 
The  objective  function  is  defined  by  subroutine  OBJFUN,  and  the  nonlinear  constraints  are  defined 
by  subroutine  CONFUN.  On  every  call,  these  subroutines  must  return  appropriate  values  of  the 
objective  and  nonlinear  constraints  in  OBJF  and  C.  The  user  should  also  provide  the  available 
partial  derivatives.  Any  unspecified  derivatives  are  approximated  by  finite  differences;  see  Section 
5.2  for  a  discussion  of  the  optional  parameter  Derivative  Level.  Just  before  either  OBJFUN  or 
CONFUN  is  called,  each  element  of  the  current  gradient  array  OBJGRD  or  CJAC  is  initialized  to  a 
special  value.  On  exit,  any  element  that  retains  the  given  v.alue  is  estimated  by  finite  differences. 

For  maximum  reliability,  it  is  preferable  for  the  user  to  provide  all  partial  derivatives  (see 
Chapter  8  of  Gill,  Murray  and  Wright,  1981,  for  a  detailed  discussion).  If  all  gradients  cannot  be 
provided,  it  is  similarly  advisable  to  provide  as  many  as  possible.  While  developing  the  subroutines 
OBJFUN  and  CONFUN,  the  Verify  parameter  (see  Section  5.2)  should  be  used  to  check  the  calculation 
of  any  known  gradients. 

4.1.  Subroutine  OBJFUN 

This  subroutine  must  calculate  the  objective  function  F(x)  and  (optionally)  the  gradient  g(x). 
The  specification  of  OBJFUN  is 

SUBROUTINE  OBJFUN<  MODE,  N,  X,  OBJF,  OBJGRD.  NSTATE  ) 

INTEGER  MODE,  N,  NSTATE 

REAL  OBJF 

REAL  X(N),  OBJGRD(N) 

Parameters: 

MODE  (Input)  This  parameter  is  set  by  NPSOL  to  indicate  the  values  that  must  be  assigned 

during  each  call  of  OBJFUN.  MODE  will  always  have  the  value  2  if  all  components  of  the 
objective  gradient  are  specified  by  the  user,  i.e.,  if  Derivative  Level  is  either  1  or  3 
(see  Section  5.2).  If  some  gradient  elements  are  unspecified,  NPSOL  will  call  OBJFUN 
with  MODE  =  0,  1  or  2. 

If  MODE  =  2,  compute  OBJF  and  the  available  components  of  OBJGRD. 

If  MODE  =  1,  compute  all  available  components  of  OBJGRD;  OBJF  is  not  required. 

If  MODE  =  0,  only  OBJF  needs  to  be  computed;  OBJGRD  is  ignored. 

(Output)  If  for  some  reason  you  wish  to  terminate  the  solution  of  the  current  prob¬ 
lem,  set  MODE  to  a  negative  value,  e.g.,  —1. 

N  (Input)  The  number  of  variables,  i.e.,  the  dimension  of  X.  The  actual  parameter  N 

will  always  be  the  same  Fortran  variable  as  that  input  to  NPSOL,  and  must  not  be 
altered  by  OBJFUN. 

X  (Input)  An  array  of  dimension  at  least  N  containing  the  values  of  the  variables  x  for 

which  F  must  be  evaluated.  The  array  X  must  not  be  altered  by  OBJFUN. 

OBJF  (Output)  The  computed  value  of  the  objective  function  F{x). 

OBJGRD  (Output)  The  available  components  of  the  gradient  vector  g{x),  i.e.,  OBJGRD(j)  con¬ 
tains  the  partial  derivative  dF/dxj. 

NSTATE  (Input)  If  NSTATE  =  1,  NPSOL  is  calling  OBJFUN  for  the  first  time.  This  parameter 
setting  allows  the  user  to  save  computation  time  if  certain  data  must  be  read  or 
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calculated  only  once.  If  there  are  nonlinear  constraints,  the  first  call  to  CONFUN  will 
occur  before  the  first  call  to  OBJFUN. 

4.2.  Subroutine  CONFUN 

This  subroutine  must  compute  the  nonlinear  constraint  functions  c(x)  and  (optionally)  their  gradi¬ 
ents.  (A  dummy  subroutine  CONFUN  must  be  provided  if  all  constraints  are  linear.)  The  t-th  row  of 
the  Jacobian  matrix  CJAC  is  the  vector  Vcj  =  (dci/dxiydcifdx2, . . .  ,dci/dx„)^.  The  specification 
of  CONFUN  is 

SUBROUTINE  C0NFUN(  MODE.  NCNLN.  N.  NROWJ, 

NEEDC,  X,  C,  CJAC,  NSTATE  ) 

INTEGER  MODE,  NCNLN,  N,  NROWJ 

INTEGER  NEEOCC*) 

REAL  X(N),  C(*),  CJAC(NROWJ,*) 

Parameters: 

MODE  (Input)  This  parameter  is  set  by  NPSOL  to  indicate  the  values  that  must  be  assigned 

during  each  call  of  CONFUN.  MODE  will  always  have  the  value  2  if  all  elements  of  the 
Jacobian  are  available,  i.e.,  if  Derivative  Level  is  either  2  or  3  (see  Section  5.2). 
If  some  elements  of  CJAC  are  unspecified,  NPSOL  will  call  CONFUN  with  MODE  =  0,  1, 
or  2: 

If  MODE  =  2,  only  the  elements  of  C  corresponding  to  positive  values  of  NEEDC  need  to 
be  set  (and  similarly  for  the  available  components  of  the  rows  of  CJAC). 
If  MODE  =  1,  the  available  components  of  the  rows  of  CJAC  corresponding  to  positive 
values  in  NEEDC  must  be  set.  Other  rows  of  CJAC  and  the  array  C  will  be 
ignored. 

If  MODE  =  0,  the  components  of  C  corresponding  to  positive  values  in  NEEDC  must  be 
set.  Other  components  and  the  array  CJAC  are  ignored. 

I 

(Output)  If  for  some  reason  you  wish  to  terminate  the  solution  of  the  current  prob¬ 
lem,  set  MODE  to  a  negative  value,  e.g.,  —1. 

NCNLN  (Input)  The  number  of  nonlinear  constraints,  i.e.,  the  dimension  of  C.  The  acturd 

parameter  NCNLN  is  the  same  Fortran  variable  as  that  input  to  NPSOL,  and  must  not 
be  altered  by  CONFUN. 

N  (Input)  The  number  of  variables,  i.e.,  the  dimension  of  X.  The  actual  parameter  N 

is  the  same  Fortran  variable  as  that  input  to  NPSOL,  and  must  not  be  altered  by 
CONFUN. 

NROWJ  (Input)  The  leading  dimension  of  the  array  CJAC.  NROWJ  must  be  at  least  1  and  at 

least  NCNLN. 

NEEDC  (Input)  An  array  that  specifies  the  indices  of  the  elements  of  C  or  CJAC  that  must 

be  evaluated  by  CONFUN.  NEEDC  need  not  be  checked  if  the  user  always  provides  all 
values,  since  the  unneeded  values  are  ignored. 

X  (Input)  An  array  of  dimension  at  least  N  contadning  the  values  of  the  variables  X  for 

which  the  constraiints  must  be  evaluated.  X  must  not  be  altered  by  CONFUN. 

C  (Output)  An  array  of  dimension  at  least  NCNLN  that  contains  the  appropriate  values 

of  the  nonlinear  constraints.  If  NEEDC(t)  >  0  and  NODE  =  0  or  2,  the  value  of  the  t-th 
constraint  at  X  must  be  stored  in  C(t).  (The  other  components  of  C  are  ignored.) 
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CJAC  (Output)  A  real  array  of  declared  dimension  (NROWJ,*),  where  the  second  dimen¬ 

sion  must  be  at  least  N,  containing  the  appropriate  elements  of  the  Jacobian  matrix 
evaluated  at  X.  (See  the  discussion  of  MODE  and  CJAC  above.) 

The  parameter  NSTATE  has  the  same  meaning  as  for  OBJFUM. 


4.3.  Constant  Jacobian  elements 

If  all  constraint  gradients  (Jacobian  elements)  are  known  (i.e.,  Derivative  Level  =  2  or  3;  see 
Section  5.2),  any  constant  elements  may  be  assigned  to  CJAC  one  time  only  at  the  start  of  the 
optimization.  An  element  of  CJAC  that  is  not  subsequently  assigned  in  CONFUN  will  retain  its  initial 
value  throughout.  Constant  elements  may  be  loaded  into  CJAC  either  before  the  call  to  NPSOL 
or  during  the  the  first  call  to  CONFUN  (signalled  by  the  value  NSTATE  =1).  The  ability  to  preload 
constants  is  useful  when  many  Jacobian  elements  are  identically  zero,  in  which  case  CJAC  may  be 
initiabzed  to  zero  and  non-zero  elements  may  be  reset  by  CONFUN. 

Note  that  constant  nonzero  elements  do  affect  the  values  of  the  constraints.  Thus,  if  CJAC(t,  j) 
is  set  to  a  constant  value,  it  need  not  be  reset  in  subsequent  calls  to  CONFUN,  but  the  value 
CJAC(i,  j)*X(j)  must  nonetheless  be  added  to  C(*). 

It  must  be  cmpheisized  that,  if  Derivative  Level  <  2,  unassigned  elements  of  CJAC  are  not 
treated  as  constant;  they  are  estimated  by  finite  differences,  at  non-trivial  expense.  If  the  user  does 
not  supply  a  value  for  Difference  Interval  (see  Section  5.2),  an  interval  for  each  component  of  x 
is  computed  automatically  at  the  start  of  the  optimization.  The  automatic  procedure  can  usually 
identify  constant  elements  of  CJAC,  which  are  then  computed  once  only  by  finite  differences. 
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5.  OPTIONAL  INPUT  PARAMETERS 

Several  optional  parameters  in  NPSOL  define  choices  in  the  problem  specification  or  the  algorithm 
logic.  In  order  to  reduce  the  mimber  of  formal  parameters  of  NPSOL,  these  optional  parameters 
have  associated  default  values  (see  Section  5.2)  that  are  appropriate  for  most  problems.  Therefore, 
the  user  needs  to  specify  only  those  optional  parameters  whose  values  are  to  be  different  from  their 
default  values.  The  remainder  of  this  section  can  be  skipped  by  users  who  wish  to  use  the  default 
values  for  all  optional  parameters.  A  complete  list  of  optional  parameters  and  their  default  values 
is  given  in  Section  5.3. 

Each  optional  parameter  is  defined  by  a  single  character  string  of  up  to  72  characters,  including 
one  or  more  items.  The  items  associated  with  a  given  option  must  be  separated  by  spaces  or  equal 
signs  (=).  Alphabetic  characters  may  be  upper  or  lower  case.  The  string 

Print  level  »  5 

is  an  example  of  an  optional  parameter. 

For  each  option,  the  string  contains  the  following  items. 

1.  The  keyword  (required  for  all  options). 

2.  A  phrase  (one  or  two  words)  that  qualifies  the  keyword  (only  for  some  options). 

3.  A  number  that  specifies  either  an  INTEGER  or  a  REAL  value  (only  for  some  options). 

Such  numbers  may  be  up  to  16  contiguous  characters  in  Fortran  77’s  I,  F,  E  or  D 
formats,  terminated  by  a  space. 

Blank  strings  and  comments  are  ignored  and  may  be  used  to  improve  readability.  A  comment  begins 
with  an  asterisk  (*)  and  all  subsequent  characters  are  ignored.  If  the  string  is  not  a  comment  and 
is  not  recognized,  a  warning  message  is  printed  on  the  specified  output  device  (see  Section  8.5). 
Synonyms  are  recognised  for  some  of  the  keywords,  and  abbreviations  may  be  used. 

The  following  are  examples  of  valid  option  strings  for  NPSOL: 

NOLIST 
warm  start 
COLD  START 

Verify  Constraint  gradients 
Start  OBJECTIVE  check  at  variable  9 

Stop  constraint  check  at  variable  =  20  *  The  is  optional 
Linear  Feasibility  tolerance  l.OE-8  *  for  IBM  in  double  precision. 

CRASH  TOLERANCE  =  .002 

*  This  string  will  be  completely  ignored. 

Hessian  Yes 
Iteration  limit  100 


5.1.  Specifleation  of  the  optional  parameters 
Optional  parameters  may  be  specified  in  two  ways,  as  follows. 

e  Using  subroutine  NPFILE  and  an  external  file 

The  subroutine  NPFILE  provided  with  the  NPSOL  package  will  read  options  from  an  external 
options  Rle,  and  should  be  called  before  a  call  to  NPSOL.  Each  line  of  the  options  file  defines  a 
single  optional  parameter.  The  file  must  begin  with  Begin  and  end  with  End.  (An  options  file 
consisting  only  of  these  two  lines  corresponds  to  supplying  no  options.) 
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The  specification  of  NPFILE  is 

SUBROUTINE  NPFILE (  lOPTNS,  INFORM  ) 

INTEGER  lOPTNS,  INFORM 

lOPTNS  must  be  the  unit  number  of  the  options  file,  in  the  range  [0,99],  and  is  unchanged  on  exit 
from  NPFILE.  INFORM  need  not  be  set  on  entry.  On  return,  INFORM  will  be  0  if  the  file  is  a  valid 
options  file  and  lOPTNS  is  in  the  correct  range.  INFORM  will  be  set  to  1  if  lOPTNS  is  out  of  range, 
and  will  be  set  to  2  if  the  file  does  not  begin  with  Begin  or  end  with  End. 

An  example  of  a  valid  options  file  is 

Begin 

Print  level  =  B 
Verify  Objective  Gradients 
End 


The  call 


CALL  NPFILE (  5,  INFORM  ) 


will  read  an  options  file  on  unit  5. 


e  Using  subroutine  NPOPTN 

The  second  method  of  setting  the  optional  parameters  is  through  a  series  of  calls  to  the  subroutine 
NPOPTN  provided  with  the  NPSOL  package.  The  specification  of  NPOPTN  is 

SUBROUTINE  NPOPTN (  STRING  ) 

CHARACTER* (e)  STRING 

STRING  must  be  a  single  valid  option  string  (see  above),  and  will  be  unchanged  on  exit.  NPOPTN 
must  be  called  once  for  every  optional  parameter  to  be  set.  An  example  of  a  call  to  NPOPTN  is 

CALL  NPOPTN (  ’Print  level  =  5’  ) 


•  Use  of  the  Nolist  and  Defaults  option 

In  general,  each  user-specified  optional  parameter  is  printed  as  it  is  read  or  defined.  By  using  the 
special  parameter  Nolist,  the  user  may  suppress  this  printing  for  a  given  call  of  NPSOL.  To  take 
effect.  Nolist  must  be  the  first  parameter  specified  in  the  options  file;  for  example 

Begin 

Nolist 

Verify  Objective  Gradients 
End 

Alternatively,  the  first  call  to  NPOPTN,  before  or  after  a  call  to  NPSOL,  must  be 

CALL  NPOPTNf  ’Nolist’  ). 

All  parameters  not  specified  by  the  user  are  automaticaUy  set  to  their  default  values.  Any 
optional  parameters  that  are  set  by  the  user  are  not  altered  by  NPSOL,  and  hence  changes  to  the 
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options  are  cuinuiative.  For  example,  calling  NP0PTN(  'Print  level  =  B’  )  sets  the  print  level 
to  5  for  all  subsequent  calls  to  NPSOL  until  it  is  reset  by  the  user.  The  only  exception  to  this 
rule  is  permitted  by  the  special  optional  parameter  Defaults,  whose  effect  is  to  reset  all  optional 
parameters  to  their  default  values  (see  Section  5.3).  For  example,  in  the  following  situation 

CALL  NPSOL  (  ...  ) 

C 

CALL  NP0PTN(  'Print  level  5'  ) 

CALL  NP0PTN(  'Iteration  limit  =100'  ) 

CALL  NPSOL  (  . . .  ) 

C 

CALL  NP0PTN(  'Defaults'  ) 

CALL  NPSOL  (  . . .  ) 

the  first  and  last  runs  of  NPSOL  will  occur  with  the  default  parameter  settings.  However,  in  the 
second  run,  the  print  level  and  iteration  limit  are  altered. 

5.2.  Description  of  the  optional  parameters 

The  following  list  (in  alphabetical  order)  gives  the  valid  options.  For  each  option,  we  give  the 
keyword,  any  essential  optional  qualifiers,  the  default  value,  and  the  definition.  The  minimum 
valid  abbreviation  of  each  keyword  is  underlined.  If  no  characters  of  an  optional  qualifier  are 
underlined,  the  qualifer  may  be  omitted.  The  letter  a  denotes  a  phrase  (character  string)  that 
qualifies  an  option.  The  letters  i  and  r  denote  INTEGER  and  REAL  values  required  with  certain 
options.  The  number  «  is  a  generic  notation  for  machine  precision,  and  denotes  the  relative 
precision  of  the  objective  function  (the  optional  parameter  Function  Precision;  see  below). 

Central  Difference  Interval  r  Default  values  are  computed 

If  the  algorithm  switches  to  central  differences  because  the  forward-difference  approximation  is  not 
sufficiently  accurate,  the  value  of  r  is  used  as  the  difference  interval  for  every  component  of  x. 
The  use  of  finite-differences  is  discussed  further  below  under  the  optional  parameter  Difference 
Interval. 

Cold  Start  Default  =  Cold  Start 

Warm  Start 

This  option  controls  the  specification  of  the  initial  working  set  in  both  the  procedure  for  finding 
a  feasible  point  for  the  linear  constraints  and  bounds,  and  in  the  first  QP  subproblem  thereafter. 
With  a  Cold  Start,  the  first  working  set  is  chosen  by  NPSOL  based  on  the  values  of  the  variables 
and  constraints  at  the  initial  point.  Broadly  speaking,  the  initial  working  set  will  include  equality 
constraints  and  bounds  or  inequality  constraints  that  violate  or  “nearly”  satisfy  their  bounds 
(within  Crash  Tolerance;  see  below).  With  a  Warm  Start,  the  user  must  set  the  ISTATE  array 
and  define  CLAMDA  and  R  as  discussed  in  Section  3.  ISTATE  values  associated  with  bounds  and 
linear  constraints  determine  the  initial  working  set  of  the  procedure  to  find  a  feasible  point  with 
respect  to  the  bounds  and  linear  constraints.  ISTATE  values  associated  with  nonlinear  constraints 
determine  the  initial  working  set  of  the  first  QP  subproblem  after  such  a  feasible  point  has  been 
found.  NPSOL  will  override  the  user’s  specification  of  ISTATE  if  necessary,  so  that  a  poor  choice  of 
the  working  set  will  not  cause  a  fatal  error.  A  warm  start  will  be  advantageous  if  a  good  estimate  of 
the  initial  working  set  is  available — for  example,  when  NPSOL  is  called  repeatedly  to  solve  related 
problems. 
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Crash  Tolerance  r  D<‘fault  =  .01 

This  value  is  used  in  conjunction  with  the  optional  parameter  Cold  start  (the  default  value). 
When  making  a  cold  start,  the  QP  algorithm  in  NPSOL  must  select  an  initial  working  set.  When 
r  >  0,  the  initial  working  set  will  include  (if  possible)  bounds  or  general  inequality  con.^traints  that 
lie  within  r  of  their  bounds.  In  particular,  a  constraird  of  the  form  a^x  >  1  will  be  included  in  the 
initial  working  set  if  \ajx  —  1|  <  r(l  +  (H).  If  r  <  0  or  r  >  1,  the  default  vrdue  is  used. 

Derivative  Level  i  Default  =  3 

This  parameter  indicates  which  derivatives  are  provided  by  the  user  in  subroutines  OBJFUN  and 
CONFUN.  The  possible  choices  for  i  are  the  following. 

t  Meaning 

3  All  objective  and  constraint  gradients  are  provided  by  the  user. 

2  All  of  the  Jacobian  is  provided,  but  some  components  of  the  objective  gradient  are 
not  specified  by  the  user. 

1  All  elements  of  the  objective  gradient  are  known,  but  some  elements  of  the  Jacobian 
matrix  are  not  specified  by  the  user. 

0  Some  elements  of  both  the  objective  gradient  and  the  Jacobian  matrix  are  not  specified 
by  the  user. 

The  value  t  =  3  should  be  used  whenever  possible,  since  NPSOL  is  more  rebable  and  will  usually 
be  more  efficient  when  all  derivatives  are  exact. 

If  I  =  0  or  2,  NPSOL  will  estimate  the  unspecified  components  of  the  objective  gradient, 
using  finite  differences.  The  computation  of  finite-dilTerence  approximations  usually  increases  the 
total  run-time,  since  a  call  to  OBJFUN  is  required  for  each  unspecified  element.  Furthermore,  less 
accuracy  can  be  attained  in  the  solution  (see  Chapter  8  of  Gill,  Murray  and  Wright,  1981,  for  a 
discussion  of  bmiting  accuracy). 

If  i  =  0  or  1,  NPSOL  will  approximate  unspecified  elements  of  the  Jacobian.  One  call  to 
CONFUN  is  needed  for  each  variable  for  which  partial  derivatives  are  not  available.  For  example,  if 
the  Jacobian  has  the  form 

/*  *  *  *  \ 

*11* 

*  *  1  * 

\  ♦  «  4>  *  f 

where  indicates  an  element  provided  by  the  user  and  indicates  an  unspecified  element, 
NPSOL  will  call  CONFUN  twice:  once  to  estimate  the  missing  element  in  column  2,  and  again  to 
estimate  the  two  missing  elements  in  column  3.  (Since  columns  1  and  4  are  known,  they  require 
no  calls  to  CONFUN.) 

At  times,  central  differences  are  used  rather  than  forward  differences,  in  which  case  twice  as 
many  calls  to  OBJFUN  and  CONFUN  are  needed.  (The  switch  to  central  differences  is  not  under  the 
user’s  control.) 

Difference  Interval  r  Default  values  are  computed 

This  option  defines  an  interval  used  to  estimate  gradients  by  finite  differences  in  the  following 
circumstances: 

1.  For  verifying  the  objective  and/or  constrmnt  gradients  (see  the  description  of  Verify,  below). 
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2.  For  estimating  unspecified  elements  of  the  objective  gradient  or  the  Jacobian  matrix. 

In  general,  a  derivative  with  respect  to  the  j-th  variable  is  approximated  using  the  interval  Sj,  where 
6j  =  r(l  +  |ijD,  with  x  the  first  point  feasible  with  respect  to  the  bounds  and  linear  constraints.  If 
the  functions  arc  well  scaled,  the  resulting  derivative  approximation  should  be  accurate  to  0(r).  See 
Gill,  Murray  and  Wright  (1981)  for  a  discussion  of  the  accuracy  in  finite-difference  approximations. 

If  a  difference  interval  is  not  specified  by  the  user,  a  finite-difference  interval  will  be  computed 
automatically  for  each  variable  by  a  procedure  that  requires  up  to  six  calls  of  CONFUN  and  OBJFUN 
for  each  component.  This  option  is  recommended  if  the  function  is  badly  scaled  or  the  user  wishes 
to  have  NPSOL  determine  constant  elements  in  the  objective  and  constraint  gradients  (see  the 
descriptions  of  CONFUN  and  OBJFUN  in  Section  4). 

Feasibility  Tolerance  r  Default  =  ^/t 

The  scalar  r  defines  the  maximum  acceptable  absolute  violations  in  linear  and  nonlinear  constraints 
at  a  “feasible”  point;  i.e.,  a  constraint  is  considered  satisfied  if  its  violation  does  not  exceed  r. 
If  r  <  0,  the  default  value  is  used.  Using  this  keyword  sets  both  optional  parameters  Linear 
Feasibility  Tolerance  and  Nonlinear  Feasibility  Tolerzmce  to  r.  (Additional  details  are 
given  below  under  the  descriptions  of  these  parameters.) 

Function  Precision  r  Default  =  e°® 

This  parameter  defines  e„,  which  is  intended  to  be  a  measure  of  the  accuracy  with  which  the 
problem  functions  F  and  c  can  be  computed.  The  value  of  £«  should  reflect  the  relative  precision 
of  1-f  |F(i)|;  i.e.,  acts  as  a  relative  precision  when  jF|  is  large,  and  as  an  absolute  precision  when 
1F|  is  small.  For  example,  if  F(x)  is  typically  of  order  1000  and  the  first  six  significant  digits  are 
known  to  be  correct,  an  appropriate  value  for  e„  would  be  l.OE-6.  In  contrast,  if  F(x)  is  typically 
of  order  10"*  and  the  first  six  significant  digits  are  known  to  be  correct,  an  appropriate  value  for 
would  be  l.OE-10.  The  choice  of  Cr  can  be  quite  complicated  for  badly  scaled  problems;  see 
Chapter  8  of  Gill,  Murray  and  Wright  (1981)  for  a  discussion  of  scaling  techniques.  The  default 
value  is  appropriate  for  most  simple  functions  that  are  computed  with  full  accuracy.  However;; 
when  the  accuracy  of  the  computed  function  values  is  known  to  be  significantly  worse  than  full 
precision,  the  value  of  €„  should  be  large  enough  so  that  NPSOL  will  not  attempt  to  distinguish 
between  function  values  that  differ  by  less  than  the  error  inherent  in  the  calculation. 

Hessian  No  Default  =  No 

Hessian  Yes 

This  option  controls  the  contents  of  the  upper-triangular  matrix  R  (see  Section  3).  NPSOL  works 
exclusively  with  the  transformed  and  re-ordered  Hessian  Hq  (5),  and  hence  extra  computation 
is  required  to  form  the  Hessian  itself.  If  Hessian  =  No,  R  contains  the  Cholesky  factor  of  the 
transformed  and  re-ordered  Hessian.  If  Hessiiui  =  Yes,  the  Cholesky  factor  of  the  approximate 
Hessian  itself  is  formed  and  stored  in  R.  The  user  should  select  Hessian  =  Yes  if  a  warm  start 
will  be  used  for  the  next  call  to  NPSOL. 

I^inite  Bound  Size  r  Default  =  10'® 

If  r  >  0,  r  defines  the  “infinite”  bound  BIGBND  in  the  definition  of  the  problem  constraints.  Any 
upper  bound  greater  than  or  equal  to  BIGBND  will  be  regarded  as  plus  infinity  (and  similarly  for  a 
lower  bound  less  than  or  equal  to  —BIGBND).  If  r  <  0,  the  default  value  is  used. 

Infinite  Step  Size  r  Default  =  max(BIGBND,  10'®) 

If  r  >  0,  r  specifies  the  magnitude  of  the  change  in  variables  that  is  treated  as  a  step  to  an 
unbounded  solution.  If  the  change  in  x  during  an  iteration  would  exceed  the  value  of  Infinite 
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Step,  the  objective  function  is  considered  to  be  unbounded  below  in  the  feasible  region.  If  r  <  0, 
the  default  value  is  used. 

Iteration  Limit  t  Default  =  max(50, 3(n  +  m^)  +  lOm^.) 

See  Major  Iteration  Limit  below. 

Linear  Feasibility  Tolerance  ti  Default  =  y/e 

Nonlinear  Feasibility  Tolerance  r2  Default  =  y/e 

The  scalars  ri  and  define  the  maximum  acceptable  absolute  violations  in  linear  and  nonbnear 
constraints  at  a  “feasible”  point;  i.e.,  a  linear  constraint  is  considered  satisfied  if  its  violation  does 
not  exceed  rj,  and  similarly  for  a  nonbnear  constraint  and  r2.  The  default  values  are  used  if  ri  or 
r2  is  non-positive. 

On  entry  to  NPSOL,  an  iterative  procedure  is  executed  in  order  to  find  a  point  that  satisfies  the 
bnear  constraint  and  bounds  on  the  variables  to  within  the  tolerance  rj.  All  subsequent  iterates 
will  satisfy  the  bnear  constraints  to  within  the  same  tolerance  (unless  ri  is  comparable  to  the 
finite-difference  interval). 

For  nonbnear  constraints,  the  feasibibty  tolerance  r2  defines  the  largest  constraint  violation 
that  is  acceptable  at  an  optimal  point.  Since  nonbnear  constraints  are  generaUy  not  satisfied 
until  the  final  iterate,  the  value  of  Nonlinear  Feasibility  Tolerance  acts  as  a  partial  termi¬ 
nation  criterion  for  the  iterative  sequence  generated  by  NPSOL  (see  the  discussion  of  Optimality 
Tolerance). 

These  tolerances  should  reflect  the  precision  of  the  corresponding  constraints.  For  example, 
if  the  variables  and  the  coefficients  in  the  linear  constraints  are  of  order  unity,  and  the  latter  are 
correct  to  about  6  decimal  digits,  it  would  be  appropriate  to  specify  rj  as  10“*. 

Linesearch  Toleranca  r  Default  =  0.9 

The  value  r  (0  <  r  <  1)  controls  the  accuracy  with  which  the  step  a  taken  during  each  iteration 
approximates  a  minimum  of  the  merit  function  along  the  search  direction  (the  smaUer  the  value 
of  r,  the  more  accurate  the  bnesearch).  The  default  value  r  =  0.9  requests  an  inaccurate  search, 
and  is  appropriate  for  most  problems,  particularly  those  with  any  nonbnear  constraints. 

If  there  are  no  nonbnear  constraints,  a  more  accurate  search  may  be  appropriate  when  it  is 
desirable  to  reduce  the  number  of  major  iterations — for  example,  if  the  objective  function  is  cheap 
to  evaluate,  or  if  a  substantial  number  of  gradients  are  unspecified. 

Major  Iteration  Limit  i  Default  =  max(50, 3(n  -b  m^)  -b  lOm^) 

Iteration  Limit 

Iters 

Itns 

The  value  of  t  specifies  the  maximum  number  of  major  iterations  allowed  before  termination. 
Setting  i  =  0  and  Major  Print  Level  >  0  means  that  the  workspace  needed  wiU  be  computed 
and  printed,  but  no  iterations  will  be  performed. 

Major  Print  Level  i  Default  =  10 

Print  Level 

The  value  of  t  controls  the  amount  of  printout  produced  by  the  major  iterations  of  NPSOL.  (See 
also  Minor  Print  Level,  below).  The  levels  of  printing  available  are  indicated  below. 
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i  Output 

0  No  output. 

1  The  final  solution  only. 

5  One  line  of  output  for  each  majoi.-  iteration  (no  printout  of  the  final  solution). 

>  10  The  final  solution  an<l  one  line  of  output  for  each  iteration. 

>  20  At  each  major  iteration,  the  objective  function,  the  Euclidean  norm  of  the 

nonlinear  constraint  violations,  the  values  of  the  nonlinear  constraints  (the 
array  c),  the  values  of  the  linear  constraints  (the  array  Aj^x),  and  the  current 
values  of  the  variables  (the  array  i). 

>  30  At  each  major  iteration,  the  diagonal  elements  of  the  matrix  T  associated  with 

the  TQ  factorization  (4)  of  the  QP  working  set,  and  the  diagonal  elements  of 
R,  the  triangular  factor  of  the  transformed  and  re-ordered  Hessian  (5). 

Minor  Iteration  Limit  i  Default  =  max(50, 3(n -I- m^, -t- m^v)) 

The  value  of  i  specifies  the  maximum  number  of  iterations  for  the  optimality  phase  of  each  QP 
subproblem. 

Minor  Print  Level  t  Default  =  0 

The  value  of  i  controls  the  amount  of  printout  produced  by  the  minor  iterations  of  NPSOL,  i.e.,  the 
iterations  of  the  quadratic  programming  algorithm.  (See  also  Major  Print  Level,  above.)  The 
following  levels  of  printing  are  available. 

t  Output 

0  No  output. 

1  The  final  QP  solution. 

5  One  line  of  output  for  each  minor  iteration  (no  printout  of  the  final  QP  solu¬ 

tion). 

>  10  The  final  QP  solution  and  one  brief  line  of  output  for  each  minor  iteration. 

>  20  At  each  minor  iteration,  the  current  estimates  of  the  QP  multipliers,  the  current 

estimate  of  the  QP  search  direction,  the  QP  constraint  values,  and  the  status 
of  each  QP  constraint. 

>  30  At  each  minor  iteration,  the  diagonal  elements  of  the  matrix  T  associated  with 

the  TQ  factorization  (4)  of  the  QP  working  set,  and  the  diagonal  elements  of 
the  Cholesky  factor  R  of  the  transformed  Hessian  (5). 

Nonlinear  Feasibility  Tolerance  r  Default  = 

See  Linear  Feasibility  Toleremce, above. 

Optimality  Tolerance  r  Default  = 

The  parameter  r  (e^  <  »■  <  1)  specifies  the  accuracy  to  which  the  user  wishes  the  final  iterate  to 
approximate  a  solution  of  the  problem.  Broadly  speaking,  r  indicates  the  number  of  correct  figures 
desired  in  the  objective  function  at  the  solution.  For  example,  if  r  is  10  ®  and  NPSOL  terminates 
successfully,  the  final  value  of  F  should  have  approximately  six  correct  figures. 
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NPSOL  will  terminate  successfully  if  the  iterative  sequence  of  i-values  is  judged  to  have  con¬ 
verged  and  the  final  point  satisfies  the  first-order  Kuhn-Tucker  conditions  (see  Section  2).  The 
serpieuce  of  iterates  is  considered  to  have  converged  at  x  if 


«11p!I  <  +  INII),  (15a) 

where  p  is  the  search  direction  and  a  the  step  length  from  (2).  An  iterate  is  considered  to  satisfy 
the  first-order  conditions  for  a  minimum  if 


<  \/^(l  +  max(l  -1-  |F(x)|,  ||gFB||))  (156) 

and 

|reSj  i</<o/  for  all  ji,  (15c) 

where  ^  is  the  projected  gradient  (see  Section  2),  ppR  is  the  gradient  of  F{x)  with  respect 
t(j  the  free  variables,  resj  is  the  violation  of  the  j-th  active  nonlinear  constraint,  and  ftol  is  the 
Nonlinear  Feasibility  Tolerance. 


Objective  Check  At  Variable  k 
Start  Constraint  Check  At  Variable  k 
St^p  Objective  Check  At  Variable  I 
Stop  Constraint  Check  At  Variable  I 


Default  =  1 
Default  =  1 
Default  =  n 
Default  =  n 


These  keywords  take  effect  only  if  Verify  level  >  0  (see  below).  They  may  be  used  to  control 
the  verification  of  gradient  elements  computed  by  subroutines  OBJFUN  and  CONFUN.  For  example, 
if  the  first  30  components  of  the  objective  gradient  appeared  to  be  correct  in  an  earlier  run,  so 
that  only  component  31  remains  questionable,  it  is  reasonable  to  specify  Start  Objective  Check 
At  Column  31.  If  the  first  30  variables  appear  linearly  in  the  objective,  so  that  the  corresponding 
gradient  elements  are  constant,  the  above  choice  would  also  be  appropriate. 


Verify  Level  t 

Verify  No 

Verify  Level  -1 

Verify  Level  0 

Verify  Objective  Gradients 

Verify  Level  1 


Default  =  0 


Verify  Constraint  Gradients 

Verify  Level  2 

Verify 

Vbrify  Yes 

Verify  Gradients 

\^rify  Level  3 

Th  esc  keywords  refer  to  finite-difference  checks  on  the  gradient  elements  computed  by  the  user- 
provided  subroutine.s  OBJFUN  and  CONFUN.  (Unspecified  gradient  components  are  not  checked.)  It 
is  possible  to  specify  Verify  Levels  0-3  in  several  ways,  as  indicated  above.  For  example,  the 
nonlinear  objective  gradient  (if  any)  will  be  verified  if  either  Verify  Objective  or  Verify  Level 
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1  is  specified.  Similarly,  the  objective  and  the  constraint  gradients  will  be  verified  if  Verify  Yes 
or  Verify  Level  3  or  Verify  is  specified. 

IfO  <  i  <  3,  gradients  will  be  verified  at  the  first  point  that  satisfies  the  linear  constraints  and 
bounds.  If  i  —  0,  only  a  “cheap”  test  will  be  performed,  requiring  one  call  to  OBJFUN  and  one  call 
to  CONFUN.  If  1  <  t  <  3,  a  more  reliable  (but  more  expensive)  check  will  be  made  on  individual 
gradient  components,  within  the  ranges  specified  by  the  Start  and  Stop  keywords  described  above. 
A  result  of  the  form  “OK”  or  “BAD?”  is  printed  by  NPSOL  to  indicate  whether  or  not  each  component 
appears  to  be  correct. 

If  10  <  t  <  13,  the  action  is  the  same  as  for  t  —  10,  except  that  it  will  take  place  at  the 
user-specified  initial  value  of  x. 

We  suggest  that  Verify  Level  3  be  specified  whenever  a  new  function  routine  is  being  de¬ 
veloped. 

5.3.  Optional  parameter  checklist  and  default  values 

For  easy  reference,  the  following  sample  NPOPTN  list  shows  all  valid  keywords  and  their  default 
values.  The  default  options  Frmetion  Precision,  Linear  Feasibility  Tolerance,  Nonlinear 
Feasibility  Tolerance  and  Optimality  Tolerance  depend  upon  e,  the  relative  precision  of  the 
machine  being  used.  The  values  given  here  correspond  to  double  precision  arithmetic  on  IBM 
360  and  370  systems  and  their  successors  (c  «  2.22  x  10“'®).  Similar  values  would  apply  to  any 
machine  having  about  16  decimal  digits  of  precision. 


*  List  of  optional  parameters. 

Central  Difference  Interval  ? 

Cold  Start 

Crash  Tolerance  .01 

Derivative  Level  3 

Difference  Interval  ? 

Function  Precision  8.2E-15 

Hessian  No 

Infinite  Boimd  l.OE+10 

Infinite  Step  1.0E-*-10 

Linear  Feasibility  Tolerance  1.5E-8 

Linesearch  Tolerance  0.9 

Major  Iteration  Limit  60 

Major  Print  Level  10 

Minor  Iteration  Limit  60 

Minor  Print  Level  0 

Nonlinear  Feasibility  Tolerance  1.5E-8 
Optimality  Tolerance  6.4E-12 

Start  Objective  Check  1 

Start  Constraint  Check  1 

Stop  Objective  Check  ? 

Stop  Constraint  Check  ? 

Verify  Level  0 


♦  Computed  automatically 

* 

* 

* 

*  Computed  automatically 

* 
e 

*  Plus  infinity 
e 

* 

* 

♦  or  3(ti -f  mt )  +  lOm^v 

e 

♦  or  3(n  -f-  mj.  4-  m^) 

♦ 

*  \fi 

*  eO.8 

e 

* 

*  n 

*  n 

*  Cheap  test 
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6.  DESCRIPTION  OF  THE  PRINTED  OUTPUT 

The  level  of  printed  output  from  NPSOL  is  controlled  Ly  the  user  (see  the  descriptiiuis  of  Hajor 
Print  Level  and  Minor  Print  Level  in  Section  5.2).  If  Minor  Print  Level  >  II,  output  is 
obtained  from  the  subroutines  that  solve  the  QP  subproblem.  For  a  detailed  description  of  this 
information  the  reader  should  refer  to  the  user’s  guide  for  LSSOL  (Gill  ct  al.,  1986a). 

When  Major  Print  Level  >  5,  the  following  line  of  output  is  produced  at  every  major 
iteration  of  NPSOL.  In  all  cases,  the  values  of  the  quantities  printed  are  those  in  effect  on  completion 
of  the  given  iteration. 

Itn  is  the  iteration  count. 

is  the  sum  of  the  iterations  required  by  the  feasibility  and  optimality  phases 
of  the  QP  subproblem.  Generally,  ItQP  will  be  1  in  the  later  iterations,  since 
theoretical  analysis  predicts  that  the  correct  active,  set  will  be  identified  near 
the  solution  (see  Section  2). 

Note  that  ItQP  may  be  greater  than  the  Minor  Iteration  Limit  if  some  it¬ 
erations  are  required  for  the  feasibility  phase. 

is  the  step  a  taken  along  the  computed  search  direction.  On  reasonably  well- 
behaved  problems,  the  unit  step  will  be  taken  as  the  solution  is  approached. 

is  the  cumulative  number  of  evaluations  of  the  objective  function  needed  for 
the  linesearch.  Evaluations  needed  for  the  estimation  of  the  gradients  by  finite 
differences  are  not  included.  Nfun  is  printed  as  a  guide  to  the  amount  of  work 
required  for  the  linesearch. 

is  the  value  of  the  augmented  Lagrangian  merit  function  (11)  at  the  current 
iterate.  This  function  will  decrease  at  each  iteration  unless  it  was  necessary 
to  increase  the  penalty  parameters  (see  Section  2.2).  As  the  solution  is  ap¬ 
proached,  Merit  will  converge  to  the  value  of  the  objective  function  at  the 
solution. 

If  the  QP  subproblem  does  not  have  a  feasible  point  (signified  by  “I”  at  the 
end  of  the  current  output  line),  the  merit  function  is  a  large  multiple  of  the 
constraint  violations,  weighted  by  the  penalty  parameters.  During  a  sequence 
of  major  iterations  with  infeasible  subproblems,  the  sequence  of  Merit  values 
will  decrease  monotonically  until  either  a  feasible  subproblem  is  obtained  or 
NPSOL  terminates  with  INFORM  =  3  (no  feasible  point  could  be  found  for  the 
nonlinear  constraints). 

If  no  nonlinear  constrsdnts  are  present  (i.e.,  NCNLN  =  0),  this  entry  contains 
Objective,  the  value  of  the  objective  function  F(x).  The  objective  function 
will  decrease  monotonically  to  its  optimal  value  when  there  are  no  nonlinear 
constraints. 

is  the  number  of  simple  bound  constraints  in  the  predicted  active  set. 
is  the  number  of  general  linear  constraints  in  the  predicted  active  set. 

is  the  number  of  nonlinear  constraints  in  the  predicted  active  set  (not  printed 
if  NCNLN  is  zero). 

is  the  number  of  columns  of  Z  (see  Section  2.1).  The  value  of  Nz  is  the  number 
of  variables  minus  the  number  of  constraints  in  the  predicted  active  set;  i.e., 
Nz  =  N  -  (End  -I-  Lin  +  Nln). 


End 

Lin 

Nln 

Nz 


ItQP 

Step 

Nfun 

Merit 
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'V'' 

i 


V, 


I 


I 


N.’. 

^l‘' 


I 


Norm  Gf 

Norm  Gz 

Cond  H 
Cond  Hz 

Cond  T 

Norm  C 

Penalty 

Conv 


M 

I 

C 


is  the  Euclidean  norm  of  the  gradient  of  the  objective  function  with  respect 
to  the  free  variables,  i.e.,  variables  not  currently  held  at  a  bound. 

is  Euclidean  norm  of  the  projected  gradient  (see  Section  2.1). 

Norm  Gz  will  be  approximately  zero  in  the  neighborhood  of  a  solution. 

is  a  lower  bound  on  the  condition  number  of  the  Hessian  approximation  H. 

is  a  lower  bound  on  the  condition  number  of  the  projected  Hessian  approxima¬ 
tion  Hz  {Hz  =  Z^Hf^Z  =  R^Rz\  see  (5)  and  (10)  in  Section  2).  The  larger 
this  number,  the  more  difficult  the  problem. 

is  a  lower  bound  on  the  condition  number  of  the  matrix  of  predicted  active 
constraints. 

is  the  Euclidean  norm  of  the  residuals  of  constraints  that  are  violated  or  in  the 
predicted  active  set  (not  printed  if  NCNLN  is  zero).  Norm  C  will  be  approximately 
zero  in  the  neighborhood  of  a  solution. 

is  the  Eucbdean  norm  of  the  vector  of  penalty  parameters  used  in  the  aug¬ 
mented  Lagrangian  merit  function  (not  printed  if  NCNLN  is  zero). 

is  a  three-letter  indication  of  the  status  of  the  three  convergence  tests  (15a)- 
(15c)  defined  in  the  description  of  the  optional  parameter  Optimality  Toler¬ 
ance  in  Section  5.  Each  letter  is  “T”  if  the  test  is  satisfied,  and  “F”  otherwise. 
The  three  tests  indicate  whether:  (a)  the  sequence  of  iterates  has  converged; 
(b)  the  projected  gradient  (Norm  Gz)  is  sufficiently  small;  and  (c)  the  norm  of 
the  residuals  of  constraints  in  the  predicted  active  set  (Norm  C)  is  small  enough. 

If  any  of  these  indicators  is  “F”  when  NPSOL  terminates  with  INFORM  =  0,  the 
user  should  check  the  solution  carefully. 

is  printed  if  the  quasi-Newton  update  was  modified  to  ensure  that  the  Hessian 
approximation  is  positive-definite  (see  Section  2.3). 

is  printed  if  the  QP  subproblem  has  no  feasible  point. 

is  printed  if  central  differences  were  used  to  compute  the  unspecified  objective 
and  constraint  gradients.  If  the  value  of  Step  is  zero,  the  switch  to  central 
differences  was  made  because  no  lower  point  could  be  found  in  the  bnesearch. 
(In  this  case,  the  QP  subproblem  is  re-solved  with  the  central-difference  gra¬ 
dient  and  Jacobian.)  If  the  value  of  Step  is  non-zero,  central  differences  were 
computed  because  Norm  Gz  and  Norm  C  imply  that  X  is  close  to  a  Kuhn- Tucker 
point. 


When  Major  Print  Level  =  1  or  Major  Print  Level  >  10,  the  summary  printout  at  the 
end  of  execution  of  NPSOL  includes  a  listing  of  the  status  of  every  variable  and  constraint.  Note 
that  default  names  are  assigned  to  aJl  variables  and  constraints. 

The  following  describes  the  printout  for  each  variable. 

Variable  gives  the  name  (VARBL)  and  index  j  {j  =  1  to  N)  of  the  variable. 

State  gives  the  state  of  the  variable  in  the  predicted  active  set  (FR  if  neither  bound  is 

in  the  active  set,  EQ  if  a  fixed  variable,  LL  if  on  its  lower  bound,  UL  if  on  its  upper 
bound).  If  the  variable  is  predicted  to  lie  outside  its  upper  or  lower  bound  by 
more  than  the  feasibility  tolerance.  State  will  be  “++”  or  “ — ”  respectively. 
(The  latter  situation  can  occur  only  when  there  is  no  feasible  point  for  the 
bounds  and  linear  constraints.) 
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The  printout  for  general  constraints  is  the  same  as  for  variables,  except  for  the  following: 
Linear  constr  is  the  name  (LNCON)  and  index  t  (t  =  1  to  NCLIN)  of  a  linear  constraint. 

Nonlnr  constr  is  the  name  (NLCON)  and  index  t  (t  =  1  to  NCNLN)  of  a  nonlinear  constraint. 


Upper  boimd 


Value  is  the  value  of  the  variable  at  the  final  iteration. 

Lower  bound  is  the  lower  bound  specified  for  the  variable.  (“None”  indicates  that  BL(j)  < 
-BIGBND.) 

Upper  boimd  is  the  upper  bound  specified  for  the  variable.  (“None”  indicates  that  BU(j)  > 
BIGBND.) 

Lagr  multiplier  is  the  value  of  the  Lagrange  multiplier  for  the  associated  bound  constraint.  This 
will  be  zero  if  State  is  FR.  If  X  is  optimal,  the  multiplier  should  be  non-negative 
if  State  is  LL,  and  non-positive  if  State  is  UL. 

Residual  is  the  difference  between  the  variable  “Value”  and  the  nearer  of  its  bounds 

BL(j)  and  BU(j). 
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7.  INTERPRETATION  OF  THE  RESULTS 

Tlie  input  data  for  NPSOL  should  always  be  checked  (even  if  NPSOL  terminates  with  the  value 
INFORM  =0!).  Two  common  sources  of  error  are  uninitialized  variables  and  incorrect  gradients, 
which  may  cause  underflow  or  overflow  on  some  machines.  The  user  should  check  that  all  compo¬ 
nents  of  A,  BL,  BU  and  X  arc  defined  on  entry  to  NPSOL,  and  that  OBJFUN  and  CONFUN  compute  all 
relevant  components  of  OBJGRD,  C  and  CJAC. 

In  the  following,  we  bst  the  different  ways  in  which  NPSOL  is  terminated  and  discuss  what 
further  action  may  be  necessary. 

Termination  Discussion  and  Recommended  Action 

Underflow  A  single  underflow  will  always  occur  if  machine  constants  are  computed  automat¬ 
ically  (as  in  the  distributed  version  of  NPSOL;  see  Section  8).  Other  floating-point 
underflows  may  occur  occasionally,  but  can  usually  be  ignored. 

Overflow  If  the  printed  output  before  the  overflow  error  contains  a  warning  about  seri¬ 

ous  ill-conditioning  in  the  working  set  when  adding  the  j-th  constraint,  it  may 
be  possible  to  avoid  the  difficulty  by  increasing  the  magnitude  of  the  optional 
parameter  Linear  Feasibility  Toleremce  or  Nonlinear  Feasibility  Toler- 
emce,  and  rerunning  the  program.  If  the  message  recurs  even  after  this  change,  the 
offending  linearly  dependent  constraint  (with  index  “j”)  must  be  removed  from 
the  problem.  If  overflow  occurs  in  one  of  the  user-supplied  routines  (e.g.,  if  the 
nonlinear  functions  involve  exponentials  or  singularities),  it  may  help  to  specify 
tighter  bounds  for  some  of  the  variables  (i.e.,  reduce  the  gap  between  appropriate 
tj  and  Uj).  If  overflow  continues  to  occur  for  no  apparent  reason,  contact  the 
authors  at  Stanford  University. 

INFORM  =  0  The  iterates  have  converged  to  a  point  X  that  satisfies  the  first-order  Kuhn-Tucker 
conditions  to  the  accuracy  requested  by  the  optional  parameter  Optimality  tol¬ 
erance  (see  Section  5.2),  i.e.,  the  projected  gradient  and  active  constraint  residuals 
are  negligible  at  X. 

The  user  should  check  whether  the  following  four  conditions  are  satisfied;  (i)  the 
final  value  of  Norm  Gz  is  significantly  less  than  that  at  the  starting  point;  (ii) 
during  the  final  major  iterations,  the  values  of  Step  and  ItQP  are  both  one;  (Hi) 
the  last  few  values  of  both  Norm  Gz  and  Norm  C  become  small  at  a  fast  bnear  rate; 
and  (iv)  Cond  Hz  is  small.  If  all  these  conditions  hold,  X  is  almost  certainly  a  local 
minimum  of  NP.  (See  Section  9  for  a  specific  example.) 

INFORM  =  1  The  point  X  satisfies  the  Kuhn-Tucker  conditions  to  the  accuracy  requested,  but 
the  sequence  of  iterates  has  not  yet  converged.  NPSOL  was  terminated  because  no 
further  improvement  could  be  made  in  the  merit  function. 

This  value  of  INFORM  may  occur  in  several  circumstances.  The  most  common 
situation  is  that  the  user  asks  for  a  solution  with  accuracy  that  is  not  attainable 
with  the  given  precision  of  the  problem  (as  specified  by  Function  Precision;  see 
Section  5.2).  This  condition  will  also  occur  if,  by  chance,  an  iterate  is  an  “exact” 
Kuhn-Tucker  point,  but  the  change  in  the  variables  was  significant  at  the  previous 
iteration.  (This  situation  often  happens  when  minimizing  very  simple  functions, 
such  as  quadratics.) 

If  the  four  conditions  listed  above  for  INFORM  =  0  are  satisfied,  X  is  likely  to  be  a 
solution  of  NP  regardless  of  the  value  of  INFORM. 


INFORM  =  2 


INFORM  =  3 


INFORM  =  4 


INFORM  =  6 


NPSOL  has  terminated  without  finding  a  feasible  point  for  tlu'  bnear  constraints 
and  bounds,  which  means  that  no  feasible  jjoint  exists  for  the  given  value  of  Linear 
Feasibility  Toleraince.  The  user  should  check  that  there  are  no  constraint 
redundancies.  If  the  data  for  the  constraints  are  accurate  only  to  an  absolute 
precision  o',  the  user  should  ensure  that  the  value  of  the  optional  parameter  I.inear 
Fea'^ibility  Tolerance  is  greater  than  a.  For  example,  if  all  elements  of  A  are  of 
oruer  unity  and  are  accurate  to  only  three  decimal  places.  Linear  Feasibility 
Tolerance  should  be  at  least  10“ 

There  has  been  a  sequence  of  QP  subproblems  for  which  no  feasible  point  could 
be  found  (indicated  by  “I”  at  the  end  of  each  terse  line  of  output).  This  behavior 
will  occur  if  there  is  no  feasible  point  for  the  nonlinear  constraints.  (However, 
there  is  no  general  test  that  can  determine  whether  a  feasible  point  exists  for  a  set 
of  nonlinear  constraints.)  If  the  infeasible  subproblems  occur  from  the  very  first 
major  iteration,  it  is  highly  likely  that  no  feasible  point  exists.  If  infeasibilities 
occur  when  earlier  subproblems  have  been  feasible,  small  constraint  inconsistencies 
may  be  present.  The  user  should  check  the  vabdity  of  constraints  with  negative 
values  of  ISTATE.  If  the  user  is  convinced  that  a  feasible  point  does  exist,  NPSOL 
should  be  restarted  at  a  different  starting  point. 

If  the  algorithm  appears  to  be  making  progress,  Major  Iteration  Limit  may  be 
too  small.  If  so,  increase  its  value  and  rerun  NPSOL  (possibly  using  the  Warm 
Start  option).  If  the  algorithm  seems  to  be  “bogged  down”,  the  user  should  check 
for  incorrect  gradients  or  ill-conditioning  as  described  below  under  INFORM  =  6. 

Note  that  ill-conditioning  in  the  working  set  is  sometimes  resolved  automatically 
by  the  algorithm,  in  which  case  performing  additional  iterations  may  be  helpful. 
However,  ill-conditioning  in  the  Hessian  approximation  tends  to  persist  once  it 
has  begun,  so  that  allowing  additional  iterations  without  altering  R  is  usually  in¬ 
advisable.  If  the  quasi-Newton  update  of  the  Hessian  approximation  was  modified 
during  the  latter  iterations  (i.c,,  an  “M”  occurs  at  the  end  of  eacii  terse  line),  it 
may  be  worthwhile  to  try  a  warm  start  at  the  final  point  as  suggested  above. 

A  sufficient  decrease  in  the  merit  function  could  not  be  attained  during  the  final 
Unesearch.  This  sometimes  occurs  because  an  overly  stringent  accuracy  has  been 
requested,  i.e..  Optimality  Tolerance  is  too  small.  In  this  case  the  user  should 
apply  the  four  tests  described  under  INFORM  =  0  above  to  determine  whether 
or  not  the  final  solution  is  acceptable  (see  Gill,  Murray  and  Wright,  1981,  for  a 
discussion  of  the  attainable  accuracy). 

If  many  iterations  have  occurred  in  which  e.sscntially  no  progress  has  been  made, 
or  NPSOL  has  failed  completely  to  move  from  the  initial  point,  subroutines  OBJFUN 
or  CONFUN  may  be  incorrect.  The  user  should  refer  to  the  comments  below  under 
INFORM  =  7  and  check  the  gradients  using  the  Verify  parameter.  Unfortunately, 
there  may  be  small  errors  in  the  objective  and  constraint  gradients  that  cannot 
be  detected  by  the  verification  process.  Finite-difference  approximations  to  first 
derivatives  are  catastrophically  affected  by  even  small  inaccuracies.  An  indication 
of  this  situation  is  a  dramatic  alteration  in  the  iterates  if  the  finite-difference 
interval  is  altered.  One  might  also  suspect  this  type  of  error  if  a  switch  is  made  to 
central  differences  even  when  Norm  Gz  and  Norm  C  arc  large. 

Another  possibility  is  that  the  search  direction  has  become  inaccurate  because  of 
ill-conditioning  in  the  Hessian  approximation  or  the  matrix  of  constraints  in  the 
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working  set;  either  form  of  ill-conditioning  tends  to  be  reflected  in  large  values  of 
ItQP  (the  number  of  iterations  required  to  solve  each  QP  subproblem). 

If  the  condition  estimate  of  the  projected  Hessian  (Cond  Hz)  is  extremely  large, 
it  may  be  worthwhile  to  rerun  NPSOL  from  the  final  point  with  the  Warm  Start 
option.  In  this  situation,  ISTATE  should  be  left  unaltered  and  R  should  be  reset  to 
the  identity  matrix. 

If  the  matrix  of  constraints  in  the  working  set  is  ill-conditioned  (i.e.,  Cond  T  is 
extremely  large),  it  may  be  helpful  to  run  NPSOL  with  a  relaxed  value  of  the 
Feasibility  Tolerance.  (Constraint  dependencies  are  often  indicated  by  wide 
variations  in  size  in  the  diagonal  elements  of  the  matrix  T,  whose  diagonals  will 
be  printed  for  Major  Print  Level  >  30.) 

INFORM  =  7  Large  errors  were  found  in  the  derivatives  of  the  objective  function  and/or  nonlin¬ 
ear  constraints.  This  value  of  INFORM  will  occur  if  the  verification  process  indicated 
that  at  least  one  gradient  or  Jacobian  component  had  no  correct  figures.  The  user 
should  refer  to  the  printed  output  to  determine  which  elements  are  suspected  to 
be  in  error. 

As  a  first  step,  the  user  should  check  that  the  code  for  the  objective  and  constraint 
values  is  correct — for  example,  by  computing  the  function  at  a  point  where  the 
correct  value  is  known.  However,  care  should  be  taken  that  the  chosen  point  fully 
tests  the  evaluation  of  the  function.  It  is  remarkable  how  often  the  values  i  =  0  or 
X  =  1  are  used  to  test  function  evalnation  procedures,  and  how  often  the  special 
properties  of  these  numbers  make  the  test  meaningless. 

Special  care  should  be  used  in  this  test  if  computation  of  the  objective  function 
involves  subsidiary  data  communicated  in  COMMON  storage.  Although  the  first 
evaluation  of  the  function  may  be  correct,  subsequent  calculations  may  be  in  error 
because  some  of  the  subsidiary  data  has  accidentally  been  overwritten. 

Errors  in  programming  the  function  may  be  quite  subtle  in  that  the  function 
value  is  “almost”  correct.  For  example,  the  function  may  not  be  accurate  to  full 
precision  because  of  the  inaccurate  calculation  of  a  subsidiary  quantity,  or  the 
limited  accuracy  of  data  upon  which  the  function  depends.  A  common  error  on 
machines  where  numerical  calculations  are  usually  performed  in  double  precision 
is  to  include  even  one  single-precision  constant  in  the  calculation  of  the  function; 
since  some  compilers  do  not  convert  such  constants  to  double  precision,  half  the 
correct  figures  may  be  lost  by  such  a  seemingly  trivial  error. 

INFORM  =  9  An  input  parameter  is  invalid.  The  user  should  refer  to  the  printed  output  to 
determine  which  parameter  must  be  re-defined. 
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8.  IMPLEMENTATION  INFORMATION 

8.1.  Format  of  the  distribution  tape 

The  source  code  and  example  program  for  NPSOL  are  distributed  on  a  magnetic  tape  containing 
12  files.  The  tape  characteristics  are  described  in  a  document  accompanying  the  tape;  normally 
they  are  9  track,  ICOO  bpi,  unlabeled,  ASCII,  80-character  records  (card  images),  4800-character 
blocks. 

The  following  is  a  list  of  the  files  and  a  summary  of  their  contents.  For  reference  purposes  we 
give  a  name  to  each  file.  However,  the  names  will  not  be  recorded  on  unlabeled  tapes.  The  MACH, 
LSCODE  and  NPCODE  files  are  composed  of  several  smaller  source  files  described  in  Section  8.3. 


File 

Name 

Type 

Cardsf 

Description 

1. 

DPMACH 

FORTRAN 

450 

Double-precision  source  file  1:  MCSUBS 

2. 

DPLSCODE 

FORTRAN 

8250 

Double-precision  source  files  2-5:  BLAS, . . . 

,  OPSUBS 

3. 

DPNPCODE 

FORTRAN 

6880 

Double-precision  source  files  6-8:  CHSUBS, . 

. . , SRSUBS 

4. 

DPLSMAIN 

FORTRAN 

260 

Double-precision  source  file  LSMAIN 

5. 

DPNPMAIN 

FORTRAN 

500 

Double-precision  source  file  NPMAIN 

6. 

LSMAIN 

DATA 

6 

Options  file  for  LSMAIN 

7. 

NPMAIN 

DATA 

14 

Options  file  for  NPMAIN 

8. 

SPMACH 

FORTRAN 

450 

Single-precision  source  file  1 

9. 

SPLSCODE 

FORTRAN 

8250 

Single-precision  source  files  2-5 

10. 

SPNPCODE 

FORTRAN 

6880 

Single-precision  source  files  6-8 

11. 

SPLSMAIN 

FORTRAN 

260 

Single-precision  version  of  file  4 

12. 

SPNPMAIN 

FORTRAN 

500 

Single-precision  version  of  file  5 

t  Approximate  figure. 


One  MACH,  one  LSCODE  and  one  NPCODE  file  should  be  selected  for  any  given  installation. 
DPMACH,  DPLSCODE  and  DPNPCODE  are  intended  for  machines  that  generally  require  double  precision 
computation.  Examples  include  IBM  Systems  360,  370,  3033,  3081,  etc.;  Amdahl  470,  Facom, 
Fujitsu,  Hitachi,  and  other  systems  analogous  to  IBM;  DEC  VAX  systems;  Data  General  MV/8000; 
ICL  2900  series;  recent  PRIME  systems;  DEC  Systems  10  and  20;  Honeywell  systems;  and  the 
Univac  1100  series. 

SPMACH,  SPLSCODE  and  SPNPCODE  are  intended  for  machines  for  which  single  precision  is  suit¬ 
ably  accurate  for  numerical  computation.  Examples  include  the  Burroughs  6700  and  7700  series: 
the  CDC  6000  and  7000  series  and  their  Cyber  counterparts;  and  the  Cray-1  and  Cray-2. 


8.2.  Installation  procedure 

1.  Obtain  the  appropriate  MACH,  LSCODE  and  NPCODE  files  from  the  tape. 

2.  If  necessary,  edit  the  subroutine  MCHPAR  according  to  Section  8.5. 

3.  Decide  whether  or  not  to  split  the  LSCODE  and  NPCODE  tape  files  into  source  files  BLAS  through 
SRSUBS  as  suggested  in  Section  8.3. 

4.  Compile  all  the  routines  that  were  originally  in  the  LSCODE  and  NPCODE  files  together  with 
those  from  MACH.  Run  them  in  conjunction  with  the  main  program  NPMAIN  from  either  file  5 
or  file  12.  Check  the  output  against  that  in  Section  9. 
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8.3.  Source  files 

NPSOL  has  been  written  in  ANSI  (1977)  Fortran  and  tested  on  an  IBM  3081K  computer  using  the 
IBM  Fortran  77  compiler  VS  Fortran.  Certain  unavoidable  machine  dependencies  arc  confined  to 
the  routine  MCHPAR. 

The  sourc.*  code  is  divided  into  8  logical  parts.  For  ease  of  handling,  these  are  combined  into 
the  MACH,  LSCODE  and  NPCOOE  files  on  the  distribution  tape,  but  for  subsequent  maintenance  we 
recommend  that  8  separate  files  be  kept.  In  the  description  below  we  suggest  a  name  for  each 
file  and  summarize  its  purpose.  We  then  list  the  names  of  the  Fortran  subroutines  and  functions 
involved.  The  naming  convention  should  minimize  the  risk  of  a  clash  with  user-written  routines. 

File  1.  HCSUBS  Computation  of  machine-dependent  constants.  ^ 

MCHPAR  MCEPS  MCENVl  MCENV2  HCSTOR 

File  2.  BLAS  Basic  Linear  Algebra  Subprograms  (a  subset). 

DASUM  DAXPY  DCOPY  DDOT  DNRH2  DSVAP  DSCAL  IDAMAX 

These  routines  are  functionally  similar  to  members  of  the  BLAS  package  (Lawson  et  al., 
1979).  If  possible  they  should  be  replaced  by  authentic  BLAS  routines.  (Versions  may 
exist  that  have  been  tuned  to  your  particular  machine.) 

DGEMV  DGERl 

These  routines  are  functionally  similar  to  members  of  the  Level  2  BLAS  packages  (Don- 
garra  et  al.,  1985). 

DCOND  DDIV  DDSCL  DLOAD  DNORM  DSSQ  DSWAP  ICOPY 

IDRANK  ILOAD 

These  arc  additional  utility  routines  that  could  be  tuned  to  your  machine.  DLOAD  is  used 
the  most  frequently,  to  load  a  vector  with  a  constant  value. 

DR0T3  DR0T3G  DGEAPQ  DGEQR  DGEQRP  DGRFG 

These  linear  algebra  routines  are  used  to  compute  and  update  various  matrix  factoriza¬ 
tions  in  NPSOL. 


File  3.  CMSUBS  General  utility  routines. 

CMALF  CMALFl  CMCHK  CMFEAS  CHPRT  CMQMUL  CMRSOL  CMRSWP 

CMRIMD  CMTSOL 

File  4.  LSSUBS  Least-squares  routines. 

LSADD  LSADDS  LSBNDS  LSCHOL  LSCORE  LSCRSH  LSDEL  LSDFLT 

LSFEAS  LSFILE  LSGETP  LSGSET  LSKEY  LSLOC  LSNOVE  LSMULS 

LSOPTN  LSPRT  LSSETX  LSSOL 

File  5.  OPSUBS  Option  string  handling  routines. 

OPFILE  OPLOOK  OPNUM  OPSCAN  OPTOKN  OPUPPR 


File  6.  CHSUBS  Derivative  checking  routines. 


CHCORE  CHFD 


CHKGRO  CHKJiC 
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File  7. 

NPSUBS 

NPCHKD 

Nonlinear 

NPCORE 

optimization  routines. 

NPCRSH  NPDFLT  NPFEAS 

JIPFILE 

NPGQ 

NPiqP 

NPKEY 

NPLOC 

NPMRT  NPOPTN  NPPRT 

NPSETX 

NPSRCH 

NPUPDT 

File  8. 

NPSOL 

SRSUBS 

SRCHQ 

Lincsearch  routines. 

SRCHC 

8.4.  Common  blocks 

Certain  Fortran  COMMON  blocks  a/e  used  in  the  NPSOL  source  code  to  communicate  between  sub¬ 
routines.  Their  names  are  listed  below. 


CMDEBG 

LSDEBG 

NPDEBG 

LSPARl 

LSPAR2 

NPPARl 

NPPAR2 

SOLICM 

S0L3CM 

S0L4CM 

S0L5CM 

S0L6CM 

SOLMCH 

SOLINP 

S0L4NP 

S0L5NP 

S0L6NP 

S0L7NP 

SOLILS 

S0L3LS 

SOLISV 

8.5.  Machine-dependent  subroutines 

The  routine  MCHPAR  in  the  MACH  file  may  require  modification  to  suit  a  particular  machine  or  a 
non-stanclard  application. 

At  the  beginning  of  NPSOL,  MCHPAR  is  called  to  assign  the  machine- dependent  constants  and 
the  staiulard  input  and  output  unit  numbers.  These  parameters  are  stored  in  the  array  VMACH(15) 
in  the  labeled  COMMON  block  SOLMCH,  and  are  defined  as  follows. 


WMACH(l)  is  NBASE,  the  base  of  floating-point  arithmetic. 
WMACH(2)  is  NDIGIT,  the  number  of  NBASE  digits  of  precision. 
WMACH(3)  is  EPS,  the  floating-point  precision. 

WMACH(4)  is  RTEPS,  the  square  root  of  EPSMCH. 

WMACH(5)  is  RMIN,  the  smallest  positive  floating-point  number. 
WMACH(6)  is  RTMIN,  the  square  root  of  RMIN. 

WMACH(7)  is  ’’**AX,  the  largest  positive  floating-point  number. 
WMACH(8)  is  RTMAX,  the  square  root  of  UMAX. 

WMACH(1())  is  NIN,  the  file  number  for  the  input  stream. 

WMACH(  11)  is  NOUT,  the  file  number  for  the  output  stream. 


V\  i:hiii  f  utiii-  MCHPAR.  the  machine  constants  are  set  in  one  of  two  ways,  depending  upon  the 
villi’  if  ?li’  !■.cll  al  variable  HDWIRE.  which  is  set  in-line. 

!■  HDVIRt  I-  FALSE  (the  value  set  for  the  distributed  copy  of  MCHPAR),  the  machine  constants 
.IT'  ;i  ;  -i  .1  Ilf "iii.ii icnllv  fill  the  machine  being  used.  If  HDWIRE  is  . TRUE .,  machine  constants 

ail)';  (iM.i''  ’  a  ti,,  lll\I  3ri(l  37(1  Seni-s  are  assigned  directly  to  the  elements  ofWMACH. 

IC  I  '  -I  •  ’  :iiii'  file  method  of  assigning  the  machine  constants,  you  should  note  the  following. 

IK  I  :i  I’.if  i  t  I  he  inaihine  constants  will  always  generate  a  single  arithmetic  underflow,  and 

)  ;  li  i.r  iniate  remedial  action  may  need  to  be  taken  if  your  machine  traps  underflow. 

l!  ,  ;  A  •..K  '’  iii;i|e||ieiit  the  in  line  assignment  of  machine  constants  for  a  machine  other  than 
•  'T  I!  Me  lilM  3tiU  371)  Senes,  MCHPAR  must  be  modified  as  follows. 

1  '  ii  iio'i  'll’  111  line  assignment  of  HDWIRE  from  .FALSE,  to  .TRUE.. 
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2.  Set  the  values  of  WMACH  appropriate  for  the  machine  and  precision  being  used.  The  values  of 
NBASE,  NDIGIT,  EPSMCH,  RMIN  and  RMAX  for  several  machines  are  given  in  the  following  table, 
for  both  single  and  double  precision;  RTEPS,  RTMIN  and  RTMAX  may  be  computed  using  Fortran 
statements.  The  values  NIN  and  NOUT  depend  on  the  machine  installation. 

For  each  precision,  we  give  two  values  for  EPSMCH,  RMIN  and  RMAX.  The  first  value  is  a  For¬ 
tran  decimal  approximation  of  the  exact  quantity;  use  of  this  value  in  MCHPAR  should  cause 
no  difficulty  except  in  extreme  circumstances.  The  second  value  is  the  exact  mathematical 
representation. 


Table  of  machine-dependent  parameters 


IBM  360/370 

Single 

CDC  6000/7000 

Single 

DEC  10/20 

Single 

Univac  1100 

Single 

DEC  Vax 

Single 

NBASE 

16 

2 

2 

2 

2 

NDIGIT 

6 

48 

27 

27 

24 

EPS 

9.54E-7 

16-5 

7.11E-1S 

2-47 

7.46E-9 

2-97 

1.50E-8 

2-96 

1.20E-7 

2-99 

RMIN 

l.OE-78 

16-85 

l.OE-293 

2-975 

l,0E-38 

2-199 

l.OE-38 

2-199 

l.OE-38 

2-198 

RMAX 

l.OE+75 

16«*(1-16'®) 

l.OE+322 

2l070(i_2-4«) 

l.OE+38 

2197(1.2-97) 

l.OE+38 

2197(1-2-97) 

l.OE+38 

2197(1-2-94) 

IBM  360/370 

Double 

CDC  6000/7000 

Double 

DEC  10/20 

Double 

Univac  1100 

Double 

DEC  Vax 

Double 

NBASE 

16 

2 

2 

2 

2 

NDIGIT 

14 

96 

62 

61 

56 

EPS 

2.22D-16 

16-1’ 

2.63D-29 

2-95 

2.17D-19 

2-89 

8.68D-19 

2-60 

2.78D-17 

2-55 

RMIN 

l.OD-78 

16-85 

l.OD-293 

2-975 

l.OD-38 

2-199 

l.OD-308 

2-1095 

l.OD-38 

2-198 

RMAX 

l.OD+75 

I6«’(I-16"'^) 

l.OD+322 

2l070(i_2-98) 

l.OD+38 

2197(1-2-89) 

l.OD+307 

2l02S(i_2-61) 

l.OD+38 

2127(i_2-56) 
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9.  EXAMPLE  PROBLEM 

This  section  describes  one  version  of  tlie  so-callcd  “hexagon”  problem  (a  different  formulation  is 
given  as  Problem  108  in  Hock  and  Schittkowski,  1981).  The  problem  is  to  determine  the  hexagon 
of  maximum  area  such  that  no  two  of  its  vertices  are  more  than  one  unit  apart  (the  solution  is  not 
a  regular  hexagon).  The  corresponding  sample  main  program  and  output  from  NPSOL  are  given 
in  the  Appendix. 

All  constraint  types  are  included  (bounds,  linear,  nonlinear),  and  the  Hessian  of  the  Lagrangian 
function  is  not  positive  definite  at  the  solution.  The  problem  has  nine  variables,  non-infinite  bounds 
on  seven  of  the  variables,  four  general  linear  constraints,  and  fourteen  nonlinear  constraints. 

The  objective  function  is 


F{x)  =  -12*6  +  *1*7  -  XiXj  -  *5*8  -I-  *4*9  -|-  *3*8- 
The  bounds  on  the  variables  are 


>  0,  -1  <  13  <  1,  *5  >  0,  *6  >  0,  *7  >0,  *8  <  0,  and  *9  <  0. 


Thus, 


/»  =  (  0,  -00,  -1,  -00,  0,  0,  0,  -00,  -00)^ 

rrt 

Ufl  =  (  00,  00,  1,  00,  00,  00,  00,  0,  0)  . 

The  general  linear  constraints  are 


Hence, 


t,.  = 


X2 

-  *1  >  0, 

X3 

-  X2 

IV 

0 

*3 

-  *4  >  0, 

and 

*4  - 

-  Xj  >  0. 

/ON 

1 

0 

0 

0 

0  0 

0 

<00^ 

0 

0 

-1 

1 

0 

0 

0  0 

0 

0 

and  xic  — 

00 

,  Ac 

0 

0 

0 

1 

-1 

0 

0  0 

0 

0 

00 

\o) 

N  0 

0 

0 

1 

-1 

0  0 

0 

0/ 

<00; 

The  nonlinear  constraints  are  all  of  the  form  Ci(*)  <  1,  for  i  =  1, . . . ,  14;  hence,  all  components 
of  In  are  —00,  and  all  components  of  are  1.  The  fourteen  functions  {ci(*)}  are 


Ci(l)  =  x\  -b*!, 

C3(®)  =  (X3  -  ij*  -l-ij, 

€5(1)  =  (*i  -  *5)*  +  (*6  -  X9)’, 
C7(*)  =  (*3  -  an)*  +1?* 

C9(l)  =  (*2  -  Xs)*  +  (X7  -  X9)*, 

Cll(x)  =  (*5  -  *3)*  -f  *1, 

Ci3(*)  =  (*4  -  *5)’  +  (*9  -  Xs)*, 


C2(*)  =  (*2  -  Xi)*  -f-  (*7  -  Xs)*, 
C4(*)  =  (*i  -  *4)*  +  (*6  -  X8)*, 

C6(x)  =  *1  +  *?, 

C8(*)  =  (X4  -  *2)*  +  (*8  -  X7)*, 
Cio(*)  =  (X4  -  *3)*  +  x|, 

Ci2(x)  =  *5  +  *J, 

Ci4(*)  =  *J  +  *|. 


An  optimal  solution  (to  five  figures)  is 

I  =  (  .060947,  .59765,  1.0,  .59765,  .060947,  .34377,  .5,  -.5,  -.34377  )’’, 
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and  F(i*)  =  —1.34996.  (The  uptiinal  objective  function  is  unique,  but  is  achieved  for  otlier  values 
of  X.)  Five  nonlinear  constraints  and  one  simple  bound  are  active  at  x.  The  sample  solution  output 
is  given  later  in  this  section,  following  the  sample  main  program  and  problem  definition. 

Two  calls  are  made  to  NPSOL  in  order  to  demonstrate  some  of  its  features.  For  the  first  call, 
the  starting  point  is: 

xo  =  (  .1,  .125,  .666666,  .142857,  .111111,  .2,  .25,  -.2,  -.25  )’’. 

All  objective  and  constraint  derivatives  are  specified  in  the  user-provided  subroutines  OBJFNl  and 
CONFNl,  i.e.,  the  default  option  Derivative  Level  =  3  is  used. 

On  completion  of  the  first  call  to  NPSOL,  the  optimal  \'ariables  are  perturbed  to  produce 
the  initial  point  for  a  second  run  in  which  the  problem  functions  are  defined  by  the  subroutines 
0BJFN2  and  C0NFN2.  To  illustrate  one  of  the  finite-difference  options  in  NPSOL,  these  routines  are 
programmed  so  that  the  first  six  components  of  the  objective  gradient  and  the  constant  elements  of 
the  Jacobian  matrix  are  not  specified;  hence,  the  option  Derivative  Level  =  0  is  chosen.  During 
computation  of  the  finite-difference  intervals,  the  constant  Jacobian  elements  are  identified  and 
set,  and  NPSOL  automatically  increases  the  derivative  level  to  2. 

The  second  call  to  NPSOL  illustrates  the  use  of  the  Warm  Start  option  to  utilize  the  final 
active  set,  nonlinear  multipliers  and  approximate  Hessian  from  the  first  run.  Note  that  Hessiw 
=  Yes  was  specified  for  the  first  run  so  that  the  array  R  would  contain  the  Cholesky  factor  of  the 
approximate  Hessian  of  the  Lagrangian. 

The  two  calls  to  NPSOL  illustrate  the  alternative  methods  of  assigning  default  parameters.  For 
the  first  run,  the  parameters  are  read  from  the  options  file  NPHAIN  DATA  supplied  on  the  distribution 
tape.  In  the  second  run,  the  parameters  are  modified  using  calls  to  subroutine  NPOPTN.  (There  is 
no  special  significance  in  the  order  of  these  assignments;  an  options  file  may  just  as  easily  be  used 
to  modify  parameters  set  by  NPOPTN.) 

The  results  are  typical  of  those  obtained  from  NPSOL  when  solving  well  behaved  (non-triviall 
nonlinear  problems.  The  approximate  Hessian  and  working  set  remain  relatively  well-conditioned. 
Similarly,  the  penalty  parameters  remain  small  and  approximately  constant.  The  numerical  results 
illustrate  much  of  the  theoretically  predicted  behavior  of  a  quasi-Newton  SQP  method.  As  x 
approaches  the  solution,  only  one  minor  iteration  is  performed  per  major  iteration,  and  the  “Norm 
Gz”  and  “Norm  C”  columns  exhibit  the  fast  linear  convergence  rate  mentioned  in  Sections  6  and  7. 
Note  that  the  constraint  violations  converge  earlier  than  the  projected  gradient.  The  final  values  of 
the  projected  gradient  norm  and  constraint  norm  reflect  the  limiting  accuracy  of  the  two  quantities. 
It  is  possible  to  achieve  almost  full  precision  in  the  constraint  norm  but  only  half  precision  in  the 
projected  gradient  norm.  Note  that  the  final  accuracy  in  the  nonbnear  constraints  is  considerably 
better  than  the  feasibility  tolerance,  because  the  constraint  violations  are  being  refined  during  the 
last  few  iterations  while  the  algorithm  is  working  to  reduce  the  projected  gradient  norm.  In  this 
problem,  the  constraint  values  and  Lagrange  multipliers  at  the  solution  are  “well  balanced”,  i.e., 
all  the  multipliers  are  approximately  the  same  order  of  magnitude.  This  behavior  is  typical  of  a 
well-scaled  problem. 
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1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 
43 

49 

50 

51 

52 

53 


»  FILE  NPMAIN  FORTRAN 

»  Sample  program  for  NPSOL  Version  4.0  February  1986. 

M4  444444  4  4444444  4444444444444444444444  444  444444444444  4  44  4444444444444444 


IMPLICIT  DOUBLE  PRECISIOIM A-K,0-Z) 


*  Set  fhe  declared  array  dimensions. 

*  NHOMA  =  fhe  declared  row  dimension  of  A. 

<  NR0t4J  =  fhe  declared  row  dimension  of  CJAC. 

»  NRCWH  =  the  declared  row  dimension  of  R. 

*  MAXN  =  maximum  no.  of  variables  allowed  for. 

*  MAXBND  =  maximum  no.  of  variables  ♦  linear  C  nonlinear  constrnts. 

*  Llh'ORK  =  the  length  of  the  integer  work  array. 

»  LWORK  =  tlie  length  of  the  double  precision  work  array. 


PARAMETER 

$ 

$ 


(NROMA  =  5,  NROUJ  =  20,  NROUR  =  10, 
MAXN  =  9,  LIUORK  =  70,  LWORK  =  1000, 
MAXBND  =  MAXN  4  nROUA  4  NROUJ) 


INTEGER 

INTEGER 

DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
EXTERNAL 


ISTATEi  MAXBND) 

IU0RX( LIWORK) 

A( KRONA. MAXN) 

BL(  MAXBND),  BUI  MAXBND) 

C( NROUJ),  CJAC( NROUJ. MAXN),  CLAMDAI MAXBND) 
03JGRD(MAXN),  RiNROWR.MAXN),  X(MAXN) 

WORK  (LWORK) 

0BJFN1,  0BJFN2,  CONFN1 ,  CONFN2 


PARAMETER  (ZERO  =  0.0,  ONE  =  1.0) 


a 

a 


Set  the  actual  problem  dimensions. 

N  =  the  number  of  variables. 

NCLIN  =  tha  number  of  general  linear  constraints  (may  be  0). 

NCNLN  =  the  number  of  nonlinear  constraints  (may  be  0). 

N  =9 

NCLIN  =  4 
NCNLN  =  14 

NSND  =  N  4  NCLIN  +  NCNLN 


«  - 

»  Arsign  file  numbers  and  the  data  arrays. 

*  NOUT  =  the  unit  nunber  for  printing. 

*  lOPTNS  =  the  unit  number  for  reading  the  options  file. 

«  Bounds  ,ge.  BIGBND  will  be  treated  as  plus  infinity. 
Bounds  .le.  -  BIGBN'D  will  bo  treated  as  minus  infinity. 

*  A  =  the  linear  constraint  matrix. 

*  BL  =  the  lower  bounds  on  x,  a'x  and  c(x). 

*  BU  =  the  upper  bounds  on  x,  a'x  and  c(x), 

*  X  =  the  initial  estimate  of  the  solution. 

»  - - - - 


54  KOUT  =  6 

55  lOPTNS  =  5 
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56  BIGBND  =  1.00*15 

57 

58  »  Set  the  natrix  A. 

59 

60  DO  40  J  =  1>  N 

61  DO  30  I  =  1,  NCLIN 

62  A(I.J)  =  ZERO 

63  30  C0HTINUE 

64  40  CONTINUE 

65  A( 1  .1  )  =  -ONE 

66  A( 1 >2)  =  ONE 

67  A(2,2)  =  -ONE 

68  A(2,3)  =  ONE 

69  A(3>3)  =  ONE 

70  A(3.4)  =  -ONE 

71  A(4>4)  =  ONE 

72  A(4>5)  =  -ONE 

73 

74  *  Set  the  bounds. 

75 

76  DO  50  J  =  1.  NBND 

77  BL(J)  =  -BIGB^ffl 

73  BU(J)  =  BIGBfJD 

79  50  CONTINUE 

60  BL( 1  )  =  ZERO 

61  BL(3)  =  -ONE 

62  BL(5)  =  ZERO 

83  BL(6)  =  ZERO 

84  BL(7)  =  ZERO 

65 

66  BU(3)  =  ONE 

67  BUI  8)  :  ZERO 

83  BUI  91  =  ZERO 

89 

90  «  Set  loNer  bounds  of  zero  for  all  four  linear  constraints. 

91 

92  DO  60  J  =  N+1 ,  N+NCUM 

93  BU(J)  =  ZERO 

94  60  CONTINUE 

95 

96  «  Set  (4>per  bounds  of  one  for  all  14  nonlinear  constraints. 

97 

93  DO  70  J  =  N  ♦  NCLIN  ♦  1,  NBND 

99  BUIJ)  =  ONE 

100  70  CONTINUE 

101 

102  *  Set  the  initial  estimate  of  X. 

103 

104  XU)  =  .1 

105  XI 2)  =  .125 

106  XI 3)  =  .666666 

107  XI4)  =  .142857 

108  XI5)  =  .111111 

109  XI6)  =  .2 

no  XI7)  =  .25 
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111 

t12 

t13 

114 

115  * 

116  * 
117  « 
IIS 

119 

120 
121 
122 

123 

124 

125  » 

126  * 
127  » 
123 

129 

130 

131 

132 

133 

134 

135 

136 

137 


X(8)  =  -.2 

X(9)  =  -.25 


Read  the  options  tile. 


CALL  NPFILEt  lOPTNS,  INFORM  ) 
IF  (INFORM  .NE.  0)  THEN 

WRITE  (NOUT,  3000)  INFORM 
STOP 
END  IF 


Solve  the  problem. 


CALL  NPSOL  (  N.  NCLIN,  NCNLN,  NROUA,  NROUJ.  NROUR> 
$  A.  BL.  BU. 

$  CONFNI.  OBJFN1. 

i  INFORM,  ITER.  ISTATE, 

$  C.  CJAC.  CLAKDA.  OBJF.  OBJGRD,  R.  X. 

$  IWORK.  LIWORK,  WORK,  LWORK  ) 

IF  (INFORM  .GT.  01  GO  TO  900 


t33  »  - 

139  »  The  tolloMing  is  tor  illustrative  purposes  only. 

140  a  A  second  run  solves  the  same  problem,  but  detines  the  objective 

141  »  and  constraints  via  the  subroutines  OBJFN2  and  CONFN2.  Some 

142  a  objective  derivatives  and  the  constant  Jacobian  elements  are  not 

143  a  supplied. 

144  a  Uc  do  a  warm  start  using 

145  a  ISTATE  (the  working  set) 

146  a  CLAMDA  (the  Lagrange  multipliers) 

147  a  R  (the  Hessian  approximation) 

143  a  trom  the  previous  run,  but  with  a  slightly  perturbed  starting 

149  a  point.  The  previous  option  file  must  have  specified 

150  a  Hessian  Yes 

151  a  for  R  to  bo  a  useful  approximation. 

152  a  - - - 

153 

154  DO  100  J  =  1 ,  N 

155  X(J)  =  X(J)  +  0.01 

156  100  CONTINUE 

157 

153  a  The  previous  parameters  are  retained  and  updated. 

159 


160 

CALL  NPOPTN(  ' 

Derivative  level 

O' ) 

161 

CALL  NFOPTIK  ' 

Verify 

No'  ) 

162 

CALL  KPOPTNi  ■ 

Warm  Start ' ) 

163 

CALL  NPOPTNI  ■ 

Major  iterations 

20'  ) 

164 

165 

CALL  NPOPTNI  • 

Major  print  level 

10' ) 

40 
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166 

167 

168 

169 

170 

171 

172 

173 

174 

175 

176 

177 

178 

179 
ISO 
131 
182 

183 

184 
135 
186 
187 
183 

189 

190 

191 

192 

193 

194 

195 

196 

197 

198 
(99 


CALL  NPSOL  (  N.  NCLIN.  NCNLN.  KROUA.  NROt4J>  NRONR. 
$  A>  BL>  BU> 

$  CONFN2.  OBJFN2. 

«  INFORM 1  ITER,  1ST ATE, 

%  C,  CJAC,  CUMDA,  OBJF,  OBJ6RD,  R,  X, 

i  lUORK,  LIUORK,  WORK,  LWORK  ) 

IF  (INFORM  .GT.  0)  GO  TO  900 
STOP 


Error  exit. 


900  WRITE  (NOUT,  3010)  INFORM 
STOP 

3000  FCRMATC/  '  NPFILE  terminated  Hith  INFORM  =*,  13) 

3010  FORMAT!/  '  NFSOL  terminated  Mith  INFORM  13) 

«  End  of  the  example  program  for  NPSOL. 

END 

^♦♦♦♦++++4+++*++**++++*+*44+*+*+4++**+******+**4**4*++*+*+**++*****+*+** 

SUBROUTINE  OBJFNK  MODE.  N,  X,  OBJF,  (%J6RD.  NSTATE  ) 

IMPLICIT  DOUBLE  PRECISION! A-M,0-Z) 

DOUBLE  PRECISION  X(N),  OBJ&RD(N) 


OBJFNl  computes  the  value  and  first  derivatives  of  the  nonlinear 
objective  function. 


200 

OBJF  s  -  X(2)»X(6) 

♦  X(1  )»X(7) 

?01 

4* 

X 

X 

♦  X(3)»X(8) 

2C2 

203 

OBJGRDd  )  =  X(7) 

204 

OBJGRO(2)  =  -  X(6) 

205 

0BJCn0(3)  =  -  X(7)  ♦ 

X(8) 

206 

OBJGRO(4)  =  X(9) 

207 

OBJGRD(5)  =  -  X(8) 

208 

OBJGRDCO)  =  -  X(2) 

209 

0BJGRD(7)  =  -  X(3)  ♦ 

X(1 ) 

210 

CBJGr!D(8)  =  -  X(5)  ♦ 

X(3) 

211 

CBJGRD(9)  =  X(4) 

212 

213 

RETURN 

214 

215  » 

End  of  OBJFNl. 

216 

217 

END 

X(5)»X(8) 


>(8  »44444+4+444444444444444444444444444444444444444444444444444444444444444 


219 

220  SUBROUTINE  CONFNK  MODE,  NCNLN,  N,  WIOWJ, 


APPENDIX.  SAMPLE  PROGRAM  AND  OUTPUT 


221 

222 

223 

224 

225 

226 
227  *- 

$ 

IMPLICIT 

INTEGER 

DOUBLE  PRECISION 

NEEDC,  X>  C>  CJAC>  NSTATE  1 

DOUBLE  PRECISION(A-H.O-Z) 

NEEDCt*) 

X(N),  C(»)>  CJAC(NROWJ.e) 

220  * 

CONEN1  computes  the  values  and  first  derivatives  of  the  nonlinear 

229  * 

constraints. 

230  • 

231  » 

The  zero  elements 

of  Jacobian  matrix  are  set  only  once.  This 

232  * 

occurs  ckirlng 

the 

first  call  to  C0NFN1 

(NSTATE  =11. 

»- 

234 

PARAMETER 

(ZERO  =  0.0.  TWO  =  2 

.01 

235 

236 

IF  INSTATE  .EQ.  11  THEN 

237 

238  * 

First  call 

to  CONFN1.  Set  all  Jacobian  eleocnts  to  zero. 

239  * 

N.B.  This 

Mil] 

i  only  Hork  Mith  'Derivative  Level  -  S'. 

240 

241 

DO  120  J  = 

!•  N 

242 

DO  110 

I  =  1 

i>  NCNLN 

243 

CJAC(IiJ)  s  ZERO 

244 

1 1 0  CONTINUE 

245 

120  CONTINUE 

246 

247 

END  IF 

243 

249 

IF  (NEEOCd  1 

•  ST. 

01  THEN 

250 

C(1  ) 

s 

Xdl»i*2  ♦  X(61»s2 

251 

CJACd.l  1 

s 

TWO*X( 1 1 

252 

CJAC( 1 >6  1 

= 

TWO»X(61 

253 

END  IF 

254 

255 

IF  <NEEDC(2) 

.GT. 

01  THEN 

256 

C(2) 

= 

(X(21  -  Xd  11»«2  ♦ 

(X(7)  -  X(6lliHi2 

257 

CJAC(2,1 ) 

=  - 

TMOt*(X(21  -  Xdll 

253 

CJAC(2.2) 

= 

TKO»(X(21  -  Xd  11 

259 

CJAC(2,6) 

TKO^IXITl  -  X(611 

260 

CJAC(2,7) 

= 

TWCH<(X(71  -  X(611 

261 

END  IF 

262 

263 

IF  (NEEDC(3I 

•  GT. 

01  THEN 

264 

C(3) 

= 

(X(31  -  Xd  11»*2  ♦ 

X(61i»2 

265 

CJAC(3,1  1 

=  - 

•ns'0»t(X(31  -  Xd  11 

266 

CJACt3,3) 

= 

TWO»(X(31  -  Xd  1) 

267 

CJACI3,5I 

= 

TWO*X( 6 1 

263 

END  IF 

269 

270 

IF  (NEEDC(4) 

.GT. 

01  THEN 

271 

C(4) 

s 

(Xd  1  -  X(41}»«2  ♦ 

(X(61  -  X(S}1»«2 

272 

CJAC(4,1  1 

s 

TWO»(Xd  1  -  Xt411 

273 

CJAC(4,4) 

=  - 

TK0»(Xd  1  -  Xt411 

274 

CJAC(4,6) 

z 

Th-O^IXtOl  -  X(811 

275 

CJAC(4,8) 

=  - 

TW0»(X(61  -  X(811 
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IF  CMEEDCCS)  .ST.  0)  THEN 


C(5I 

CJACtS.I ) 
CJAC(5,5) 
CJAC(5>61 
CJAC(5,9) 
END  IF 


(X(1)  -  X(5)»»«2  ♦  <X(6)  -  X(9))»«2 

TKO«(Xm  -  X<5)1 

-  TWO*(X(  1  )  -  X(5n 
TU0i>(X(6l  -  X(9n 

-  TW0»(X(6)  -  X(9n 


IF  (MEEDC(6J  .ST.  0)  THEN 
C(6)  =  X(2)««2 

CJAC(6>2)  =  THO!<X(2) 

CJAC(6.7)  =  TU’0i*X(7l 

END  IF 

IF  (NEE0C(7)  .ST.  0»  THEN 


X(7J»»2 


C(7J  = 

CJAC(7,2I  = 
CJACt7,3)  = 
CJAC(7,7)  = 

END  IF 


(X(3I  -  X(2))**2  ♦  X(7)**2 

-  TNOtKXtS)  -  XC2n 
TV!0)f(X(3)  -  X(2n 
TMe*X<7) 


IF  (NEEDC(a)  .ST.  0)  THEN 


c(ai 

CJAC(8,2) 
CJAC(a,4l 
CJAC(8,7) 
CJAC(8,8) 
END  IF 


<XC4)  -  X(2))»»2  ♦  (X<81  -  X(7))»*2 

-  TMO*(X(A)  -  X(2n 
TKO«(X(41  -  X<2)l 

-  TH0»(X(81  -  X(7n 
TH0»(X(8J  -  X<7H 


IF  (NEE0C(9J  .ST.  0)  THEN 


C(9) 

CJAC<9,2) 
CJAC(9,5) 
CJAC(9,7) 
CJAC(9,9) 
END  IF 


(X(2I  -  X(5n»«2  ♦  (X(7J  -  X(9))»«2 

TMO*(X(2>  -  X(5)) 

«  TW0«(X(2)  -  X(5n 
TMO»(X(7)  -  X(9)l 
-  TMO»(X(7l  -  X(9)l 


IF  (NEEDC(tO)  .ST.  0)  THEN 


C(»0)  = 

CJAC(10.3)  = 
CJAC(10,4)  = 
CJAC( 10.31  = 
END  IF 


(X(4)  -  X(3)J*it2  ♦  X(8)««2 

-  TWO»(X(4)  -  X(3)l 
TWO»(X(4l  -  XC3)) 

TKO»X( 8 » 


IF  (NEEOClin  .ST.  0)  THEN 

C(11  )  =  (X(5I  -  X(3)»«»2 

CJACU1,3I  =  -  TM0»(X(5I  -  X13n 
CJACm.S)  =  TWO»(X(5)  -  XC3n 

CJAC(1I,9)  =  TMO«X?9) 

END  IF 


X19)»*2 


IF  (NEE0C(12)  .ST.  0)  THEN 
C(12I  =  X(4)«»2  ♦ 


X(8»«*2 
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331 

332 

333 

334 

335 

336 

337 

338 

339 

340 

341 

342 

343 

344 

345 

346 

347 

348 

349 

350 

351  * 

352 

353 


CJAC(12,4)  =  TV(0«X(4) 

CJAC(12,6)  =  TK0«X(8) 

END  IF 


IF  1NEEDCM3) 
C(13) 

CJAC(13,4) 
CJAC(13>5) 
CJAC( 13.8) 
CJAC(13>9) 
END  IF 


.61.  0)  THEN 
=  (X(4)  -  X(5))««2  ♦ 

=  TW0»(X(4I  -  X(5)) 

=  -  TWO*(X<4)  -  Xt5)» 

=  -  TMO»(X(9)  -  X(8)) 

=  TWO»1X<9)  -  X(8)) 


(X(91 


X(BIHHi2 


IF  (NEEDC(14) 
C(  14) 

CJAC(14,5) 
CJAC(14,9) 
END  IF 


GT.  0)  THEN 
X(5)*»2  ♦ 

TNOi*X<5) 
TMO«X(  9 ) 


X19)«H»2 


RETURN 

End  of  C0NFN1. 
END 


354  *+*♦♦♦+♦♦♦♦♦♦+♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦ 

355 

356  SUBROUTINE  OBJFNEI  MODE.  N.  X.  OBJF,  OBJGRO.  NSTATE  ) 

357  IMPLICIT  DOUBLE  PRECISia«A-H,0-Z) 

35S  DOUBLE  PRECISION  X(N).  OBJGRO(N) 

359 


360  * - 

361  «  03JFN2  cotnputes  tfie  value  and  scow  first  derivatives  of  the 

362  *  nonlinear  objective  fuiction. 

363  « - 

364 

365  OBJF  =  -  X(2)»X(6)  ♦  X(1)«X(7)  -  X(3)«X(7)  -  X(5)«X(d) 


366 

367 

»  ♦  X(4)«<Xt9)  ♦  X(3)«X(a) 

368 

OBJGRD(3)  =  -  X(7) 

♦  X(8) 

369 

OBJGROI?)  =  -  X<3) 

♦  xn ) 

370 

371 

0BJGRD(8)  =  -  X(5) 

♦  X(3) 

372 

373 

RETURN 

374  * 

375 

End  of  OBJFN2. 

376 

END 

377  »♦♦♦♦♦++♦♦♦♦♦♦+♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦ 

378 

379 

SUBROUTINE  CONFN2( 

MODE.  NCNLN.  N,  NROWJ. 

300 

381 

$ 

NEEDC.  X,  C.  CJAC,  NSTATE  ) 

382 

IMPLICIT 

DOUBLE  PRECISION! A-H, 0-2) 

333 

INTEGER 

NEEOCie) 

304 

385 

DOUBLE  PRECISION 

X(N).  C(«),  CJACCNROHJ.e) 
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386  *— 

387  * 
383  » 

CONFN2  computes  the  values  and  the 
the  nonlinear  constraints. 

390 

PARAMETER 

(TMO  =  2.01 

391 

392 

IF  (NEEDCCI ) 

.GT. 

0)  THEN 

393 

C(  1  ) 

s 

X(1)«*2  ♦  X(6HH 

394 

CJACd.n 

= 

TMO»X(1  > 

395 

CJAC(1,6) 

= 

TMO*X(6) 

396 

END  IF 

397 

393 

IF  (NEEDCtZ) 

.GT. 

0)  THEN 

399 

C(2) 

= 

<X(2I  -  XM  n»«2 

400 

CJAC(2.1 I 

=  - 

TKO»(X(2»  -  X(1  n 

401 

CJAC(2>2) 

s 

Two«(X(2i  -  xm» 

402 

CJAC(2)6) 

=  - 

TKOa(X(7)  -  Xi6n 

403 

CJAC(2,7) 

= 

TKO»(X(71  -  X(6n 

404 

END  IF 

405 

406 

IF  (NEEDC(3i 

.GT. 

0)  THEN 

407 

C(3) 

= 

(Xt3)  -  xm)»*2 

403 

CJAC(3,1  » 

r  - 

TWO»(X(3)  -  X(1  n 

409 

CJAC(3,3) 

= 

TMO»(X(3»  -  xmi 

410 

CJAC(3>6I 

= 

TWO«X(61 

411 

END  IF 

412 

413 

IF  (NEE0C(4) 

.GT. 

01  THEN 

414 

C(4) 

= 

(X(1  1  -  Xt4))»*2 

415 

CJAC(4,1 I 

5 

TWO*(X(  1  1  -  X(4n 

416 

CJAC(4,4I 

=  - 

TUO»(X(  1  1  -  X(4I) 

417 

CJAC(4,6) 

s 

TMo»(X(6i  -  xten 

418 

CJAC(4,8) 

=  - 

Two*(x(6)  -  Xian 

419 

END  IF 

420 

421 

IF  (HEEOCIS) 

.GT. 

01  THEN 

422 

C(S) 

= 

(X(1  1  -  X(511»»2 

423 

CJACfS.I ) 

s 

TW(H*(X(  1  1  -  X(51) 

424 

CJAC(5>5) 

=  - 

TMO«(X(11  -  X<5II 

425 

CJAC(5,6» 

TW0»(X(61  -  X(911 

426 

CJAC(5t9) 

=  - 

TWO«(X(61  -  x«9n 

427 

END  IF 

428 

429 

IF  (NEEDC(6) 

.GT. 

01  THEN 

430 

C(6) 

s 

X(21«»2  ♦  Xt7HH 

431 

C>>AC(6>2I 

= 

TNO»X(21 

432 

CJAC(6>7) 

s 

TMO«X171 

433 

END  IF 

434 

435 

IF  (NEEDC(7) 

.GT. 

01  THEN 

436 

C(7) 

= 

(X(31  -  X(2))»«2 

437 

CJAC(7,2) 

=  - 

TWO»(X(31  -  X(21J 

438 

CJAC(7,3) 

= 

TVIO»(X(3l  -  X(2)> 

439 

CJAC(7,7) 

s 

TWO^»X(71 

440 

END  IF 

X(6>»»2 


(X<6)  -  X(d))<H(2 


X(9))»»2 


X<7)»*2 
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441 

442 

443 

444 

445 

446 

447 
443 

449 

450 

451 

452 

453 

454 

455 

456 

457 
453 

459 

460 

461 

462 

463 

464 

465 

466 

467 
46S 

469 

470 

471 

472 

473 

474 

475 

476 

477 

478 

479 

480 

481 
432 
483 
434 
465 
486 
437 
483 
4S9 

490 

491 

492 

493 

494  * 

495 


IF  (NEEDC(8)  .GT.  0)  THEN 

C(8)  =  (X(4I  -  X<2))»«2  ♦  (Xtaj  - 

CJAC(a,2)  =  -  TWO*(Xt4)  -  X(2)) 

CJAC(a,4)  =  TMO«(X(4)  -  X(2)) 

cjAC(a,7)  =  -  TMO*(X(a»  -  x(7n 

cjActa.ai  =  TMO»(X(a)  -  x(7»> 

END  IF 

IF  1NEEDC(9)  .GT.  0)  THEN 

C(9)  =  (X(2)  -  X15)>»»2  ♦  <X(7»  - 

CJAC19,2)  =  TWO*(X(2)  -  XISI) 

CJAC(9,5)  =  -  TWO»(X(2)  -  X(S)) 

CJAC(9,7)  =  TKO»(X<71  -  X(9)» 

CJAC(9,9)  =  -  TH0»t(X(7)  -  X(9)1 

END  IF 

IF  (NEEOC(IO)  .GT.  0)  THEN 

C(10)  =  (X(4)  -  X(3)I«»2  ♦  X18NW2 

CJAC(10>3I  =  -  TWO«(X(4)  -  X(3n 

CJAC(ia>4)  =  TVIO)((X(4)  -  X(3)) 

CJAC(10)8)  =  TWO*X«8) 

END  IF 

IF  (NEEDC(II)  .GT.  0)  THEN 

C(11)  =  (X(51  -  X(3n«*2  ♦  X19)**2 

CJAC(11,3)  =  -  TW0*f(X(5)  -  X13)) 

CJAC<11,5)  =  TWQ»(X(51  -  Xt3n 

CJAC(11,9)  =  TMO»X(9) 

END  IF 

IF  (NEEDC(12>  .GT.  0)  THEN 

C(12»  =  X(4)»«2  ♦  X(8)4«2 

CJAC(12,4)  =  TM0»X(4) 

CJAC(12,8)  =  TNO«X(8) 

END  IF 

IF  (NEE0C(13)  .GT.  0)  THEN 

C(13)  =  (X(4)  -  X(5)l«»2  ♦  1X19)  - 

CJAC<13,4J  =  TW0*(X14)  -  X15)) 

CJAC(13|5)  =  -  T|{0»(X(4)  -  XC5M 

CJAC(13,8)  =  -  TWO«(X(9l  -  Xtail 

CJAC(13,9)  =  TW0X(X(9»  -  X(8)) 

END  IF 

IF  (NEE0Cn4)  .GT.  0)  THEN 

C(14l  =  X(5)»»2  ♦  X19)»»2 

CJAC(14,5)  =  TW0»X(5I 

CJACI 14,9)  =  TW0XX(9) 

EIO  IF 

RETURN 


X(7))**2 


Xt9))«K2 


Xf8))<Hi2 


End  of  C0NFN2 
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APPENDIX.  SAMPLE  PROGRAM  AND  OUTPUT 


7. 


M 


47 


i 

1 .OOE-OI 

i .31E-07 

1 

2.00000048E-Ot 

2.00000040E-01 

OK 

1 .28E-07 

2 

-4.99999523E-02 

-4.99999523E-02 

OK 

i .49E-07 

3 

-1 .i3333169E*00 

-I.I3333I89E*00 

OK 

1 .38E-07 

4 

-8.571 39826E-02 

-6.571398262-02 

OK 

I.40E-07 

5 

-2.222I9229E-02 

-2.222I9229E-02 

OK 

Column 

X(  Jl 

DXiJ) 

Rom 

Jacobian  Valua 

Difference  Approxn 

2 

1.2SE-0I 

I.28E-07 

2 

4.99999523E-02 

4.99999523E-02 

OK 

i.33E-07 

6 

2.50000000E-OI 

2.50000000E-OI 

OK 

1 .49E-07 

7 

-i. 063331 94E«00 

-1. 083331 94E*00 

OK 

1 .40E-07 

8 

-3.57140303E-02 

-3.57I40303E-02 

OK 

i.43E-07 

9 

2.777602942-02 

2.77760294E-02 

OK 

lira 


to  Jacobian  oleaienio  out  of  tha  10  sat  in  cols  I  through  C  saaa  to  ba  ok. 
Tha  largest  ralattva  arror  was  2.13E-II  In  roM  9>  coluan  2 

Verification  of  tha  objactiva  gradients. 


The  objective  gradients  seen  to  ba  ok. 

Directional  derivative  of  tha  objective 
Difference  approximation 


i.20539630E-Ot 

I.20539630E-0I 


I  r  . 
I  x". 


J  XIJi 

7  2.50E-01 

8  -2.00E-01 

9  -2.50E-0I 


DXiJI 

2.26E-06 

2.17E-06 

2.26E-06 


6(J1 

-5. 0666599 7E-0i 
S.55554986E-0I 
i.428S70t5E-0l 


Olfferanca  approxn  Itns 


-5.66665997E-0I 

5.55554986E-0I 

l.428570l5E-0t 


3  Objective  gradients  out  of  the  3  sat  in  cels 
Tha  largest  relative  arror  Mas  2.21E>1I  in  aleaant 


OK 

OK 

OK 


7  through 


9  saea  to  ba  ok. 


■ 

Itn 

I  top 

step  Nfvai  Merit 

Bnd 

Lin  Nln 

Ite 

Nora  6f 

Nora  6z 

Cond  H 

Cond  Hz 

Cond  T 

Nora  C 

Penalty  Conv 

■ 

0 

5 

0.0E*00 

1  -3.i349i7E-0i 

3 

0 

i 

5 

8.0E-OI 

3. 72-01 

i.E^OO 

1.2*00 

1.2*00 

8.82-01 

0.02*00 

F  FF 

I 

9 

i.OE^OO 

2  -i.075027E«00 

1 

0 

3 

5 

2.2E«00 

l.5E*00 

1.2*02 

7.2*00 

2.E*00 

6.62-01 

1.32*00 

F  FF 

2 

4 

l.0E*00 

3  -i.26e553E«00 

i 

0 

4 

4 

i.7E«eo 

3.3E-0i 

9.2*00 

l.E*00 

2.E*00 

1.32-01 

1.32*00 

F  FF 

3 

2 

1 .0E*00 

4  -1 .33I667E*00 

1 

0 

5 

3 

i.9C*00 

E.SE-01 

4.E*0i 

2.E*00 

2.E*00 

1.12-01 

1.32*00 

F  FF 

4 

i 

l.0E«00 

5  -1 .349354E*00 

1 

0 

5 

3 

i.6E»00 

4.5E-02 

3.E*01 

1 .E*00 

2.E*00 

1.42-02 

1.32*00 

F  FF 

k  ^ 

5 

i 

1 .0E*00 

6  -1 .349874E>00 

i 

0 

5 

3 

t.6E»00 

6.72-03 

3.E*0i 

2.E*00 

2.E*00 

9.12-04 

1.32*00 

F  FF 

6 

i 

i .0E«00 

7  -i .3499I3E«00 

1 

0 

5 

3 

t.OE*00 

5.3E-03 

3.E*01 

2.2*00 

2.E*00 

5.72-05 

1.32*00 

F  FF 

7 

t 

1.0E*00 

8  -i .349963E*00 

i 

0 

5 

3 

1.8E*00 

i.2E-03 

1.2*02 

2.E*00 

2.E*00 

3.12-04 

6.82*00 

F  FF 

8 

i 

1 .0E«00 

9  -1 .349963E«00 

i 

0 

5 

3 

I.OE^OO 

1.62-04 

1.2*02 

3.E*00 

2.2*00 

9.02-07 

6.82*00 

F  FF 

9 

i 

i .0E«00 

10  -i .349963E*00 

i 

0 

5 

3 

1.6E«00 

5.4E-06 

3.E*01 

2.E*00 

2.E*00 

1.22-06 

6.82*00 

F  TT 

10 

i 

1 .0E«00 

11  -i .349963E*00 

1 

0 

5 

3 

1.0E*08 

2.0E-07 

4.2*01 

2.E*0D 

2.E*00 

6.42-11 

6.82*00 

F  TT 

1 1 

1 

t .0E*00 

12  -1 .349963E*00 

i 

0 

5 

3 

1 .8E*00 

1.  IE-00 

1.E*02 

2.E*00 

2.2*00 

4.72-14 

6.82*00 

T  TT 

*  v* 

Exit 

MP  phase.  INFORM  a  0  MAJITS  = 

11 

NFUf 

e 

12  NGRAO  s  12 

Vr 

Exit 

NPGOL 

-  Optimal 

solution  found. 

\  ‘i 


User  s  Guide  for  NPSOL  4-0 


Final  nonlinaar  objactlw  valua  = 


Calls  io  NPOPTN 


Deri  vat I vs  laval 
Verify 
Mara  Start 
Majer  iterations 
Major  print  level 


SOL/NPSOL  -  Version  4.0  Feb  i9e6 


Paraaeters 


Linear  constraints. 
Variables . 


Nonlinear  constraints.. 
Nonlinear  Jacobian  vars 
Nonliiiear  object iv  vars 


Major  iterations  liait. 
Minor  iterations  liait. 

Uorkspace  provided  is 
To  solve  problea  Ma  need 


4 

Mne»r  fcasIbllHy . 

I.49E-08 

HARM  start . 

9 

Infinite  bound  size.... 
Infinite  step  size . 

i.OOE^IO 

1.00E*i0 

Crash  tolerance . 

14 

Cptiaality  tolerance... 

5.36E-i2 

Function  precision . 

9 

Nofiftlnear  feasibt  1  i  ty.  a 

1 .49E-08 

9 

L 1 nesearch  tolerance. . . 

9.00E-0i 

2.22E-i6 

Derivative  level . 

0 

Verify  level . 

20 

Major  print  level . 

iO 

81 

Minor  print  level...... 

0 

Mi  10001. 
Hi  9681. 


11k?  user  sots  44  out  of  126  Jacobian  aleaents. 

Each  iteration?  82  Jacobian  aleaents  Mill  be  estiaated  mawrically. 

Tlie  user  sets  9  out  of  9  objective  gradient  aleaents. 

Each  iteration?  6  gradient  eleaents  Mill  be  estiaated  nuerarically. 


Coaputation  of  the  finite-difference  intervals 


XU) 

Foraard  DXU) 

Central  DXU) 

Error  est. 

7.09E-02 

1.935067E-06 

1 .935067E-05 

1.979764E-08 

6.08E-01 

2.90482iE-06 

2. 90482  IE-05 

1.3le833E-08 

1 .OOE^OO 

3.6i3750E-07 

3.613750E-06 

O.OOOOOOE^OO 

6.08E-01 

2. 904821 E-06 

2. 90482 IE-05 

1.318833E-08 

7.09E-02 

1.935067E-06 

1 .935067E-05 

1.979764E-08 

3.54E-0i 

2.446096E-06 

2.446096E-05 

I.566159E-08 

5.10>'  '’i 

2.7283812-07 

2. 728381 E-06 

O.OOOOOOE^OO 
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6  -4.90E-0! 
9  -3.3<<E-0t 


2.69ZZ44E-07 

Z.40995eE-06 


2.692Z44E-06 

2.4099SeE-05 


O.OOOOOOE*00 
1 .5&9644E-08 


02  constant  constraint  gradient  eleiBents  assigned. 
0  constant  objective  gradient  elcuents  assigned. 


All  eissing  Jacobian  eleaents  are  constants. 


Derivative  level  increased  to 


lin 

ItQP 

Step 

Nfcn 

Merit 

Bnd 

Lin 

Nln 

He 

More  6f 

More  Ez 

Cond  H 

Cond  Hz 

Cond  T 

More  C 

Penalty 

Conv 

0 

i 

0.0E*00 

1 

-1 .349i88E*00 

1 

0 

5 

3 

1.8E*00 

1 .5E-02 

3.EF01 

4.E*00 

1 .E«00 

2.eE-02 

2.2E*00 

F  FF 

1 

t 

1 .0E*00 

3 

-1 .349963E*00 

1 

0 

5 

3 

l.8E*00 

I.3E-03 

1.E*02 

7.E»00 

t .E*00 

3.0E-04 

3.0E*02 

F  FFM 

2 

1 

1 .OEFOO 

4 

-1 .349963E  +  00 

1 

0 

5 

3 

l.6E«00 

3.5E-04 

6.E*0I 

6.E«00 

2.E«00 

7.8E-07 

2.1t*01 

F  FF 

3 

1 

1 .0E*00 

5 

-1 .349963E*00 

1 

0 

5 

3 

1 .6E*00 

2.0E-04 

6.E«01 

3.E+00 

2.E+00 

2.3E-08 

7.7E*00 

F  FF 

1 

t .0E«00 

6 

-1 .349963E*00 

1 

0 

5 

3 

1.8E«00 

7.4E-06 

9.E*0i 

3.E«00 

2.E*00 

3.9E-08 

7.7E+00 

F  FF 

5 

1 

1 .OE+00 

7 

-1 .349963E«00 

1 

0 

S 

3 

1 .0E*OO 

5.9E-07 

2.E+02 

3.E«00 

2.E*00 

4.0E-1I 

7.7E+00 

F  TT 

6 

1 

t .0E«00 

0 

-1 .34996 3E*00 

1 

0 

5 

3 

1.6E«00 

2.6E-09 

6.E«0I 

2.E*00 

l.E«00 

2.0E-13 

7.7E«00 

T  TT 

Exit 

NP  phas«. 

INFORM  =  0  MAJITS  = 

6  NFUN  = 

0  NERAO  = 

7 

Variable 

State 

Value 

LoMer  bound 

U|9pcr  bound 

Lagr  eultiplier 

Residual 

VARBL 

1 

FR 

0.6094665E- 

-01 

O.OOOOOOOE^OO 

None 

O.OOOOOOOE^OO 

0.6095E-0I 

VARBL 

2 

FP 

0.5976493 

None 

None 

O.OOOOOOOE^OO 

0.1000E*I6 

VARBL 

3 

UL 

1 .000000 

-1.000000 

1 .000000 

-0.6075429 

O.OOOOE*00 

VARBL 

4 

FR 

0.5976493 

None 

None 

0.0000000E*00 

0.I000E*I6 

VARBL 

5 

FR 

0.6094665E- 

-01 

0.0000000E«00 

None 

o.oooooooE«oo 

0.6095E-OI 

VARBL 

6 

FR 

0.3437715 

O.OOOOOOCE^OO 

None 

o.oooooooE«oo 

0.3430 

VARBL 

7 

FR 

0.5000000 

O.OOOOOOOE^OO 

None 

O.OOOOOOOEFOO 

0.5000 

VARBL 

0 

FR 

-0.5000000 

None 

o.oooooooE*oo 

O.OOOOOOOEFOO 

0.5000 

VARBL 

9 

FR 

-0.3437715 

Hone 

O.OOOOOOOEFOO 

o.oooooooE«oo 

0.3430 

Linear 

'  constr 

State 

Value 

Lower  bound 

Upper  bou-id 

Lagr  ■ultipliar 

Residual 

LUCON 

t 

FR 

0.5367026 

O.OOOOOOOEFQO 

None 

O.OOOOOOOEFOO 

0.5367 

LKCON 

2 

FR 

0.4023507 

o.oooooooE«oo 

None 

O.OOOOOOOEFOO 

0.4024 

LJXON 

3 

FR 

0.4023507 

o.oooooooE*oo 

None 

O.OOOOOOOEFOO 

0.4024 

L^JCON 

4 

FR 

0.5367027 

O.OOOOOOOEFQO 

None 

O.OOOOOOOEFOO 

0.5367 

Non 1 nr 

constr 

State 

Value 

Lower  bound 

Upper  bound 

Lagr  Multiplier 

Residual 

NLCOK 

i 

FR 

0.1216933 

None 

1.000000 

O.DOOOOOOE^OO 

0.0781 

NLCON 

2 

FR 

0.3124571 

None 

1 .000000 

O.OOOOODOEFOO 

0.6675 

NLCOfJ 

3 

UL 

1  .000000 

None 

1.000000 

-0.63184068-01 

-0.1652E-12 

MLCON 

4 

UL 

1.000000 

None 

1.000000 

-0.3202625 

-0.1104E-I2 

NICON 

5 

FR 

0.4727152 

None 

1.000000 

O.OOOOOOOEFOO 

0.5273 

HLCCM4 

6 

FR 

0.6071047 

None 

1.000000 

0.00000008*00 

0.3920 

ULCOM 

7 

FR 

0.4118861 

None 

1.000000 

0.00000008*00 

0.5881 

NLCOH 

0 

UL 

1 .000000 

None 

1.000000 

-0.1992983 

0.0000E*00 

MICON 

9 

UL 

1.000000 

None 

1.000000 

-0.3202625 

-0.6882E-i4 

NLCCM 

10 

FR 

0.4118861 

None 

1.000000 

o.oooooooE*oo 

0.5601 

inccM 

11 

UL 

1.000000 

None 

1.000000 

-0.8316406E-0I 

-0.2665E-13 

MLCON 

12 

FR 

0.6071047 

None 

1 .000000 

o.oooooooE*oo 

0.3928 

tJLCOH 

13 

FR 

0.3124571 

None 

1.000000 

o.oooooooE*oo 

0.6875 

NLCOtl 

14 

FR 

0.1210933 

None 

1.000000 

0.'''1000E*00 

0.8761 

Exit  NP30L  -  Optieal  solution  found. 

Final  nonlinear  objective  value  3  -1.349963 
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INDEX 

Ai,  (gciiPral  lini'iir  ronstrciint  matrix),  1,  5,  21. 

Arw-  (.lacobiaii  of  nonliucar  coniitTaints),  2,  5  (also 
.SCO  Jacobi, -u)  matrix). 

A,  7  (definition). 

Accuracy 

desired  in  optimal  solution,  21-22,  27  (also  see 
Optimality  Tolerance), 
of  finite-dilfercnce  gradients,  19. 
of  linesearch,  20. 

of  nouliiK-.'ir  constraints  at  solution,  22,  27 
(also  see  Nonlinear  Feasibility  Tolerance), 
of  projected  gradient  at  solution,  22,  27. 
Accur.ate  linesearch,  when  appropriate,  20. 

Active  constraints 
definition,  2. 
predicted,  3,  24. 
residuals  at  solution,  8,  20,  22. 

Active  simple  bound,  2  (also  see  Fixed  variable). 
Algorithm  of  NPSOL,  description,  2-6. 
a  (step  length  in  major  iteration),  2,  4,  6,  22,  24. 
choice  of,  4,  6. 
printed  value,  24. 

Amdahl  470,  30. 

ANSI  (1977)  Fortran,  1,  30. 

Approximate 

gradients  (see  Finite-difference  approxima¬ 
tions). 

Hessian  of  Lagrangian  function,  3,  4,  0,  19,  21. 
ASCII,  30. 

Assignment  of  constant  elements  in  Jacobian,  10. 
Attainable  .accuracy,  18. 

Augmented  Lagrangian  merit  function,  4,  6,  20. 
printed  value.  24. 

Automatic  computation  of  finite-difference  inter¬ 
vals,  18  -19. 

BAD?,  23. 

B. adly  scaled  problems,  19. 

Begin  (in  options  file),  15-16. 

BFCS  (piasi-Newton  update,  4,  6  (also  see  Ap¬ 
proximate  Hessian  of  Lagrangian  function). 
BIGBHD.  8.  19  20,  26. 

BL.  7  8  (definition),  9,  26. 

BLAS.  31. 

Level  2,  31. 

Bnd,  3,  24. 

Bounds  .and  linear  constraints,  separate  treat¬ 
ment  of,  3,  4,  6,  9,  17. 

BU.  8  (definition),  9,  26. 

Burroughs  C700  and  7700,  30. 

c(z)  (nonlinear  constraints),  1,  3,  6. 

printout  of,  21. 

('  (predicted  .active  set),  2,  3. 

Cfr  ,  2,  3,  5. 

C  (array  of  nonlinear  constraints),  10  (definition), 
13. 

C  (printed  indication  of  switch  to  central  differ¬ 
ences),  25. 

CD(’  GOOD  .and  7000,'30. 

Central  Difference  Interval,  17  (definition). 
Central  differences,  switch  to,  25. 


Cheap  gradient  test,  23. 

Checklist  of  ojitional  parameters,  23. 

Cholesky  factor,  3,  5,  6,  10. 

CJAC,  10  (definition),  14. 

CLAMDA,  10  (definition),  17. 

Cold  Start,  9,  10,  17  (definition). 

Comment  (in  optional  parameter  specification), 
15. 

Common  blocks,  list  of,  32. 

Cond  H,  25. 

Cond  Hz,  25. 

Cond  T,  25. 

Conditions  for  optimality,  2,  8,  21-22,  25,  27. 
CONFUN  (user-provided  subroutine) 

calls  needed  for  unspecified  Jacobian  elements, 
18. 

definition  as  parameter  of  NPSOL,  8. 
specification,  13-14. 

Constant  Jacobian  elements 
assignment  of  10,  14. 
automatic  computation  of,  14,  35. 

Constrained  linear  least-squares,  1. 

Constrained  stationary  point  for  QP,  5. 
Constraints 

dependencies,  resolution  of,  27,  28-29. 
nonlinear,  specification  by  user  (see  CONFUN). 
status  indicator  (sec  ISTATE). 
violation,  maximum  acceptable,  19,  20  (also 
see  Linear  Feasibility  Tolerance  and  Nonlin¬ 
ear  Feasibility  Tolerance). 

Conv  (printout  of  convergence  test  status),  25. 
Convergence  test,  21-22  (also  see  Optimality 
conditions) . 

Cost 

of  automatic  computation  of  finite-difference 
intervals,  19. 

of  unspecified  objective  gradient  elements,  18. 
of  unspecified  Jacobian  elements,  18. 

Crash  Tolerance,  17,  18  (also  see  Cold  Start). 
Cray-1  and  Cray-2,  30. 

Cyber,  30. 

d  (search  direction  in  QP  method),  4-5. 

Data  General  MV/8000,  30. 

DEC  Systems  10  and  20,  30. 

DEC  VAX,  30. 

Default  values  of  optional  parameters,  checklist 
of,  23. 

Defaults  (optional  parameter),  16-17. 
Dependencies,  constraint,  resolution,  27,  28-29. 
Derivative 

checking  (see  Verify). 

finite-difference  (sec  Finite-difference  approxi¬ 
mations). 

specification  (.see  Derivative  Level). 
Derivative  Level,  3,  10,  12,  13,  14,  18,  35. 
Diagonals 

of  fZ,  printout,  21. 
of  T,  printout,  21. 

Difference  Interval,!,  14,  18. 

use  in  approximating  unspecified  gradients,  19. 
use  in  verification  of  gradients,  18. 
Discontinuities,  isolated,  1. 

Distribution  tape,  format  of,  30. 


INDEX 
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Double  ]>rccision 

table  of  machine  constants,  33. 
version  of  code,  30. 

DOUBLE,  7. 

Effects  of  ill-conditioning,  27,  28-20. 

End  (in  options  file),  15  16. 

EPS,  32. 

c  (machine  precision),  17,  10,  32,  33. 

(h  (function  precision),  17,  10,  21,  27. 

EQ  (printed  constraint  status),  25. 

Equality  constraint,  1,  8. 

Errors  in  gradients,  20. 

Estimated  Lagrange  multiplier  (see  Lagrange 
multiplier). 

Example  problem  for  NPSOL,  34-35. 

External  file,  use  for  option  specification,  15-16. 

F(x)  (objective  function),  1. 

Facom,  30. 

Failure  in  linescarch,  28. 

Feasibility  phase  in  QP  method,  4,  5,  20. 

selection  of  initial  working  set,  0. 

Feasibility  Tolerance,  10  (definition). 
Finite-difference 

approximations  to  gradients,  1,  3,  12. 
checking  of  gradients  (see  Verify), 
intervals,  automatic  computation,  14,  10. 
tradeoffs  in  computing,  18. 

First-order  Kuhn- Tucker  conditions,  2,  8,  22,  27. 
Fixed  variable,  2,  4,  8. 

Formal  parameters  of  NPSOL,  7. 

Format  of  distribution  tape,  30. 

Fortran  77,  1,  31. 

Fortran  subroutines,  naming  convention,  31-32. 
FR  (subscript),  3  (definition),  4  (also  see  IVee 
variable) . 

FK  (printed  constraints  status),  25. 

Free  variable,  2,  3,  4. 

Fujitsu,  30. 

Function  precision  (see  cr). 

Ftmctlon  Precision,  17,  10  (definition),  27. 

FX  (subscript),  3  (definition),  4. 

g{x)  (objective  gradient),  2. 

Sfr,  2. 

Gabor,  Zsa  Zsa,  10. 

Global  convergence,  6. 

Gradient 

approximations  (see  Finite-difference  approxi¬ 
mations). 

constraint  (sec  Jacobian  matrix), 
of  Lagrangian  function,  6. 
projected  (see  Projected  gradient), 
specification  by  user  (see  CONFON  and  OBJFUN). 

H  (approximate  Hessian  of  Lagrangian  function), 
3,  6,  10,  25. 

Hq,  4,  6,  10,  25. 

H,  10. 

Hz,  25. 

HDWIRE,  32. 

Hessiiui  approximation  (see  Approximate  Hessian 
of  Lagrangian  function). 


Hessian,  10  (definition). 

Hessian,  trsmsformed  and  reordered  (sc-e  Hq). 
Hexagon  example,  34. 

Hitachi,  30. 

Honeywell,  30. 

I  (printout  indicating  infeasible  QP  subproblem), 
24,  25. 

IBM 

360/370  and  3033/3081,  30. 

VS  Fortran,  31. 

ICL  2000  series,  30. 

Identity  matrix,  in  resetting  Hessian,  20. 

Dl  conditioning,  effects  of,  27,  28-20. 
Implementation  information,  30-33. 

Inaccuracies,  effect  of,  28. 

Inaccurate  lincsearch,  20. 

Inconsistent  linear  constraints,  treatment,  28. 
Incorrect  gradients,  8,  28  (also  see  Verify). 
Inequality  constraints  (nonlinear),  treatment  in 
merit  function,  6. 

Infeasible  problem 

in  QP  subproblem,  5,  0,  24,  25,  28. 

for  bounds  and  linear  constraints,  4,  8,  28. 

for  nonlinear  constraints,  8,  28. 

Infeasibilities,  4,  5,  24,  25. 

Infinite  Bound  Size  (BIGBND),  10  (definition). 
Infinite  lower  or  upper  bound,  1,  8. 

Infinite  Step  Size,  20  (definition). 

INFORM,  8  (definition). 

Initial  working  set  in  QP  subproblem, 
with  Cold  Start,  9,  17-18. 
with  Warm  Start,  0,  17. 

Input  parameter,  invalid,  8,  29. 

Installation  procedure,  30. 

Interpretation  of  results,  27-20. 

Invalid  input  parameter,  8,  29. 
lOPTNS  (options  file  number),  15-16. 

Isolated  discontinuities,  1. 

1ST  ATE,  0  (definition),  17. 

printout,  25,  26. 

ITER,  8  (definition). 

Iteration  Limit,  20  (definition). 

Iters, 20. 

Itn  (printed  value),  24. 

Itns,  20. 

ItqP  (printed  value),  24. 

IW,  11  (definition). 

Jacobian  matrix  (nonlinear  constraints),  2,  3,  8, 
10,  14. 

assignment  of  constant  elements,  10,  14. 
specification  by  user  (see  COHFUN). 
unspecified  elements,  18. 

Kuhn- Tucker  conditions,  first-order,  2,  8,  22,  27. 
Keyword  in  option  specification,  15. 

t  (lower  bound  vector),  1,  3,  7-8  (also  see  BL). 
Lack  of  progress  in  m^or  iteration,  28. 

Lagr  multiplier  (printed  valne),  26. 

Lagrange  multiplier,  2,  3,  6,  10,  26. 
of  QP  subproblem,  5,  21. 
optimal,  5. 
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A,  2,  C  (also  sec  Lagrange  multiplier). 

LENIU,  11  (definition). 

LENW,  11  (definition). 

Level  2  BLAS,  31. 

Limiting  accuracy,  18. 

Lin,  3,  24. 

Linear  constr,  26. 

Linear  Feasibility  Tolerance,  4,  0,  20  (defini¬ 
tion). 

adjustment  to  avoid  overflow,  27. 

Linear  least-squares  code  (see  LSSOL). 

Lines  of  code  in  NPSOL,  1,  30. 

Linescarcli,  4,  0,  20  (also  sec  Step  length), 
effect  of  accuracy,  20. 
routines  for,  32. 

Linesearch  Tolerance,  20  (definition). 

LL  (printed  constraint  status),  25. 

LNCON,  26. 

Local  minimum  (see  Optimality  conditions). 
Lower  bound  (in  printout),  26. 

LSSOL,  1,  3,  4. 

m  (number  of  constraints  in  predicted  active 
set),  3. 

mi.  (number  of  general  linear  constraints),  1. 
m^'  (number  of  nonlinear  constraints),  1. 

M  (printed  indicator  of  modified  Hessian  update), 
6,  25,  28. 

Machine  constants 
computation  of,  32. 
tables  of,  33. 

Machine  precision  (see  e). 

Major  iteration,  2. 

Major  Iteration  Limit,  8,  20  (definition),  28. 
Major  Print  Level,  8,  11,  20  (definition),  24,  25. 
Maximum  acceptable  constraint  violations  (see 
Linear  Feasibility  Tolerance  and  Nonlin¬ 
ear  Feasibility  Tolerance). 

MCHPAR,  32  (also  see  Machine  constants). 

Merit  function,  4,  6,  20,  24. 

Merit  (printed  value),  24. 

Method 

of  NPSOL,  description,  2-6. 

QP,  4-5. 

Minimum  abbreviation  (of  optional  parameter), 
15. 

Minimum  sum  of  infeasibilities  in  QP,  5  (also  see 
Feasibility  phase). 

Minor  iteration  (within  QP  method),  2,  3,  4-5. 
Minor  Iteration  Limit,  21  (definition). 

Minor  Print  Level,  21  (definition),  24. 

MINOS,  1. 

MODE 

in  CONFDN,  13. 
in  OBJFUN,  12. 

Modification  of  quasi-Newton  update,  6,  25,  28. 

/ijv,  6. 

Multiplier  (see  Lagrange  multiplier). 

n  (ntmibcr  of  variables),  1,  3  (also  see  N). 
npH  (number  of  free  variables),  2,  3. 
mfx  (number  of  fixed  variables),  2  (also  see  Bnd). 
nz,  3,  24. 

N,  7  (definition),  12,  13. 


Naming  convention,  Fortran  subroutines,  31-32. 
Natural  order  of  variables,  10. 

NBASE,  32. 

NCLIN,  7  (definition)  (also  see  mi). 

NCNLN,  7  (definition),  13. 

NDIGIT,  32. 

NEEOC,  13. 

Nfun  (printed  value),  24. 

NIN,  32. 

NLCON,  26. 

Nln  (printed  value),  3,  24. 

No  feasible  point 

for  bounds  and  linear  constraints,  4,  8,  28. 
for  nonlinear  constraints,  8,  28. 
in  QP  subproblcm,  5,  9,  24,  25,  28. 

No  progress  in  linesearch,  8,  28. 

Noxist  option,  16. 

Non-existent  lower  or  upper  bound,  1,  8. 

Nona  (in  printout),  26. 

NOUT,  32. 

Nonlinear  Feasibility  Tolerance,  0,  20  (defini¬ 
tion),  21. 

adjustment  to  avoid  overflow,  27. 

Nonlinear  constraints 

inequality,  in  merit  function,  6. 
predicted  active  set,  3. 
specification  by  user  (see  CONFUN). 
violated,  residuals  of,  25. 

Nonlinear  optimization,  routines  for,  32. 

Nonlnr  constr,  26. 

Norm  C,  25. 

Norm  Cf,  3,  25. 

Norm  Gz,  3,  25. 

NP  (problem  statement),  1,  2. 

HPFILE,  15-16. 

NPOPTH,  16. 

list,  sample  16. 

NPSOL 

algorithm  of,  2-6. 
lines  of  code,  1,  30. 
parameters  of,  7-11. 
specification,  7. 
solving  related  problems,  17. 

NROWA,  7  (definition). 

NROWJ,  7  (definition),  13. 

NROWS,  7  (definition). 

NSTATE,  12,  14. 

Null  space,  3. 

dimension  of  (see  nz). 

Nz,  24. 

Objective  (printed  value),  24. 

Objective  function  (f’(i)),  1. 
precision  of  (see  cn). 
specification  by  user  (see  OBJFUN). 

OBJF,  10  (definition),  12. 

OBJFUN  (user-provided  subroutine) 

calls  needed  for  unspecified  gradient  clcsments, 
18. 

definition  as  parameter  of  NPSOL,  8. 
specification,  12-13. 

OBJCRD,  10  (definition),  12. 

OK,  23. 

Optimal 


Mir 
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Lagrange  multiplier,  5. 
solution  (see  Optimality  conditions). 
Optimality 

conditions,  2,  8,  21-22,  25,  27. 
phase,  in  QP  method,  4,  5,  0,  21. 

Optimality  Tolerance,  20,  21-22  (definition),  25, 
27. 

Option-handling  routines,  31. 

Optional  parameters,  definition,  15-23. 

Options  file,  15-16. 

Ordering  of  variables,  10. 

Orthogonal  transformation,  3. 

Output  (sec  Printout). 

Overflow,  27. 

p  (search  direction  in  m^or  iteration),  2,  4. 
Parameters 

of  CONFON,  13-14. 
of  NPSOL,  7-11. 
of  OBJFUN,  12-13. 

Penalty  parameters  (in  merit  function),  6,  25. 
Penalty  (printed  value),  25. 

Phase  1  (see  Feasibility  phase). 

Phase  2  (see  Optimality  phase). 

Phrase  (to  modify  optional  parameter),  15. 
Positive-definite  Hessian  approximation  (see  Ap¬ 
proximate  Hessian  of  Lagrangian  function). 
Precision 

function  (see  e^). 
machine  (see  t). 

of  linear  constraints,  relation  to  Linear  Feasi¬ 
bility  Tolerance, 28. 

Predicted  active  set  (see  Active  constraints  and 
Working  set). 

Prcloading  constant  Jacobian  elements,  10,  14. 
Primal  method  (for  QP),  4. 

Prime  Systems,  30. 

Print  Level,  20  (definition). 

Printout 

control  of,  20-21. 
description,  24-26. 

Programming  errors,  symptoms,  20. 

Projected  gradient 

of  nonlinear  objective,  2,  3,  8,  21-22,  25. 
of  QP  subproblem,  5. 

Q,  3,  6. 

Qfr,  3,  5. 

Quadratic  program 

jiietbod  of  LSSOL,  4-5. 
multipliers,  3,  5,  6. 
subproblem,  2,  4-5. 

Qualifying  phrase  (in  optional  parruneter),  15. 
Quasi-Newton 

approximation  (sec  Approodmate  Hessian  of 
Lagrangian  function), 
update,  4,  6. 

QPSOL,  1. 

R, 3 

Rzy  3,  25. 

R,  10  (definition),  10,  21. 

Rank-one  update  to  R,  0. 


Rank-two  modification  (sec  Quasi-Newton  up¬ 
date). 

REAL,  7. 

Re-ordered  Hessian  (see  Approximate  Hessian  of 
Lagrangian  function). 

References,  36. 

Related  problems,  solved  by  NPSOL,  17. 
Resetting 

Hessian  matrix,  to  overcome  ill-conditioning, 

20. 

optional  parameters  to  defaults,  16-17. 
Residual  (printed  value),  26. 

Residuals,  constraint 

allowed  maximum  at  solution  (see  Linear  Fea¬ 
sibility  Tolerance  and  Nonlinear  Feasi¬ 
bility  Tolerance), 
in  optimality  conditions,  22 
Resolution  of  constraint  dependencies,  27,  28-20. 
Reverse-triangular  matrix,  3  (also  see  T). 
p  (see  Penalty  parameters). 

RNAX,  32. 

RMIN,  32. 

RTEPS,  32. 

RTNAX,  32. 

RTHIN,  32. 

Scaling  techniques,  10. 

Search  direction 

in  m^or  iteration,  2,  4. 
in  QP  subproblem,  4-5. 

Separate  treatment  of  bounds  and  linear  con¬ 
straints,  3,  4,  6,  0,  17. 

Sequential  quadratic  programming  algorithm  (see 
SQP  algorithm). 

<r  (step  length  in  QP  method),  5. 

Single  precision 

table  of  machine  constants,  33. 
version  of  code,  30. 

Singularities  in  objective  function,  27. 

Slack  variables  in  merit  function,  6. 

Source  files,  list,  30. 

Sparse  problems,  1. 

Specification 

of  CONFON,  13-14. 
of  NPSOL,  7. 
of  OBJFON,  12-13. 

SQP  algorithm,  2-4,  6. 

Start  Constraint  Check,  22  (definition). 

Start  Objective  Check,  22  (definition). 

State  (printed  value),  25. 

Status  of  constraints  (see  ISTATE). 

Step  (printed  value),  24. 

Step  length 

in  miyor  iteration  (a),  2,  4,  6,  22,  24. 
in  QP  method  (it),  5. 

Stop  Constraint  Check,  22  (definition). 

Stop  Objective  Check,  22  (definition). 

Sufficient  decrease  (see  Step  length). 

Sum  of  infeasibilities 
in  QP,  4,  5. 

of  nonlinear  constraints,  22,  25. 

Synonyms  (for  optional  parameters),  15. 

T,  3,  5,  21,  25. 
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Tape 

chararteristics,  30. 
format,  30. 

Termination 

criteria,  8,  20  (also  see  Optimality  conditions), 
user-controlled,  8  (see  NODE). 

TQ  factorization,  3,  5. 

TVansforined  and  re-ordercd  Hessian  (see  Ap¬ 
proximate  Hessian  of  Lagrangian  function). 
Two-phase  primid  method  for  QP,  4. 

u  (vector  of  upper  bounds),  1,  3,  8  (also  see  BO). 
OL  (printed  constraint  status),  2S. 

Unbounded  objective  function,  20. 

Underflow,  27. 

Univac  1100,  30. 

Unspecified  derivatives,  1,  18. 

Update 

of  Hessian  approximation  (see  Quasi-Newton 
update) . 

of  working  set  in  QP  method,  5. 

Updating  matrix  factorizations,  routines  for,  31. 
Upper  bound  (in  printout),  26. 

Upper-triangular  matrix  (see  Cholesky  factor). 
User-requested  termination  (see  MODE). 
User-supplied  subroutines,  12-14. 


(priutiHl  constraint  status),  25  (;ilso  see  Infea¬ 
sible  problem). 


1 


Valid  option  strings,  examples  of,  16. 

Value  (printed  value),  26. 

VIRBL,  25. 

Variable,  25. 

Verification  of  gradients,  4,  18,  22-23,  20. 
Verify,  4,  12,  22  (definition). 

Verify  Level,  22  (definition). 

Vertex,  5. 

Violations,  constraint  (see  Infeaaibilities). 

W,  11  (definition). 

Warm  start,  example  of,  35 
Warm  Start,  0,  10,  17  (definition). 

Wejl  scaled  problems,  19. 

WHICH,  32  (also  see  Machine  constants). 

Working  precision  (see  c). 

Working  set,  3,  4,  9. 
changes  in,  5. 
initial,  in  QP,  17-18. 

Condition  estimate  (see  Cond  T). 

Workspace  parameters,  11. 

X  (vector  of  unknowns),  1,  2. 
printout,  25. 

X,  11  (definition),  12,  13. 

(  (Lagrange  multipliers  for  active  bounds),  2. 

X*  (solution  of  NP),  2,  3. 

3- 

2  (basis  for  null  space),  2,  5,  24. 

Z  9fr,  2,  3. 

Zero  Jacobian  elements,  14. 

--  (printed  constraint  status),  25  (also  see  Infea¬ 
sible  problem). 
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